Assignment 4: System Planning¶

In this group work, we will investigate scenarios for a 100% renewable electricity system for the country ofSlovakia, which is a landlocked, mountainous country in Europe, lies roughly between latitudes 47° and 50° N, and longitudes 16° and 23° E. Forests cover 41% of Slovak land surface. [wikipedia].¶

In [1]:
# importing the necessary libraries & Packages:

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import geopandas as gpd
import pypsa
import cartopy.feature as cfeature
import cartopy.crs as ccrs
import cartopy
import warnings
import atlite
import rasterio
from rasterio.plot import show
import xarray as xr
import netCDF4
warnings.filterwarnings('ignore')
from atlite.gis import ExclusionContainer
from atlite.gis import shape_availability
import rasterio as rio
import os
import country_converter as coco

import logging
logging.basicConfig(level=logging.INFO)

from shapely.geometry import Polygon, MultiPolygon, Point
In [2]:
# extracting GADM data belonging to Slovakia & creating the shape of the country including its regions:

svk_regions = gpd.read_file('gadm_410-levels-ADM_1-SVK.gpkg')
svk_regions = svk_regions.set_index('NAME_1')
display(svk_regions)
shape_svk = svk_regions.to_crs(3035).geometry     # Slovokia is in Europe (crs = 3035)
shape_svk.plot()
GID_0 GID_1 COUNTRY geometry
NAME_1
Banskobystrický SVK SVK.1_1 Slovakia MULTIPOLYGON (((19.45671 48.09119, 19.45659 48...
Bratislavský SVK SVK.2_1 Slovakia MULTIPOLYGON (((17.31282 48.45557, 17.31436 48...
Košický SVK SVK.3_1 Slovakia MULTIPOLYGON (((22.10827 48.38136, 22.10702 48...
Nitriansky SVK SVK.4_1 Slovakia MULTIPOLYGON (((18.83296 47.81751, 18.81611 47...
Prešovský SVK SVK.5_1 Slovakia MULTIPOLYGON (((21.54594 48.79453, 21.54533 48...
Trenčiansky SVK SVK.6_1 Slovakia MULTIPOLYGON (((17.41709 48.83292, 17.42872 48...
Trnavský SVK SVK.7_1 Slovakia MULTIPOLYGON (((17.45999 48.14781, 17.46112 48...
Žilinský SVK SVK.8_1 Slovakia MULTIPOLYGON (((19.02072 48.87508, 19.01839 48...
Out[2]:
<Axes: >
No description has been provided for this image
In [3]:
###### !!! this section is performed only for presentation purposes !!   #########

df = pd.read_csv('global_power_plant_database.csv')     # Reading the data frame included all power plant data of all countries
df_svk = df[df.iloc[:, 0] == 'SVK']                     # extracting the Slovakia-related data
display(df_svk.head(5))
df_total_generation_groupby = df_svk.groupby('primary_fuel').estimated_generation_gwh_2017.sum()
#df_total_generation_groupby.head(5)

# Plotting the pie chart of total generation per technology in Slovakia:
plt.figure(figsize=(8, 5))
plt.pie(df_total_generation_groupby, labels=df_total_generation_groupby.index, autopct='%1.1f%%', startangle=90)
plt.legend(df_total_generation_groupby.index, bbox_to_anchor=(1, 0.5), loc="center left", title="Energy Carriers")
plt.title('Estimated Generation GWh (2017) by Primary Fuel for Slovakia')
plt.savefig('Estimated Generation(2017).png')
plt.show()
print()
print()

geometry = gpd.points_from_xy(df_svk["longitude"], df_svk["latitude"])
gdf_svk = gpd.GeoDataFrame(df_svk, geometry=geometry, crs=4326)
country country_long name gppd_idnr capacity_mw latitude longitude primary_fuel other_fuel1 other_fuel2 ... estimated_generation_gwh_2013 estimated_generation_gwh_2014 estimated_generation_gwh_2015 estimated_generation_gwh_2016 estimated_generation_gwh_2017 estimated_generation_note_2013 estimated_generation_note_2014 estimated_generation_note_2015 estimated_generation_note_2016 estimated_generation_note_2017
19861 SVK Slovakia Badín WKS0063339 6.0 48.2650 19.0560 Solar NaN NaN ... 8.29 8.91 8.90 8.58 8.36 SOLAR-V1-NO-AGE SOLAR-V1-NO-AGE SOLAR-V1-NO-AGE SOLAR-V1-NO-AGE SOLAR-V1-NO-AGE
19862 SVK Slovakia Bohunice Nuclear Power Plant Slovakia GEODB0002385 880.0 48.4960 17.6816 Nuclear NaN NaN ... NaN NaN NaN NaN 6805.78 NO-ESTIMATION NO-ESTIMATION NO-ESTIMATION NO-ESTIMATION CAPACITY-FACTOR-V1
19863 SVK Slovakia Cierny Vah Pumped Storage Hydroelectric Power ... GEODB0042635 735.0 49.0088 19.9122 Hydro NaN NaN ... 2597.78 2210.18 2927.73 2216.24 2210.18 HYDRO-V1 HYDRO-V1 HYDRO-V1 HYDRO-V1 HYDRO-V1
19864 SVK Slovakia Cunovo Hydroelectric Power Plant Slovakia GEODB0042641 24.0 48.0299 17.2254 Hydro NaN NaN ... 118.03 117.69 129.12 122.91 130.73 HYDRO-V1 HYDRO-V1 HYDRO-V1 HYDRO-V1 HYDRO-V1
19865 SVK Slovakia Dubnica nad Vahom Hydroelectric Power Plant Sl... GEODB0042662 16.5 48.9638 18.1454 Hydro NaN NaN ... 73.11 55.75 63.33 52.24 55.80 HYDRO-V1 HYDRO-V1 HYDRO-V1 HYDRO-V1 HYDRO-V1

5 rows × 36 columns

No description has been provided for this image

Out[3]:
<GeoAxes: >
No description has been provided for this image
In [4]:
# Reading the file included all of the Exclusive economic zones of the countries, which is not the case for Slovakia because it is a landlocked country:
# But the files are read anyway:

all_EEZ = gpd.read_file('eez_v11.gpkg') 
all_EEZ = gpd.read_file('eez_boundaries_v11.gpkg')
In [5]:
## !!! For presentation purposes !!! ###

# Reading the Slovakia Shape from data "country_shapes.geojson" file:

url = 'country_shapes (4).geojson'
countries = gpd.read_file(url).set_index("name")
#countries.tail(5)
gebco_svk = "GEBCO_2014_2D-SK.nc"
gebco = rasterio.open(gebco_svk)
band = gebco.read(1)
fig, ax = plt.subplots(figsize=(8, 12))
countries.loc[["SK"]].plot(ax=ax, color="none")
show(band, transform=gebco.transform, cmap="RdBu_r", ax=ax)
Out[5]:
<Axes: >
No description has been provided for this image

First, a land eligibility analysis, based on the criteria listede belo, is performed as follows:¶

Criteria for onshore wind implementation:

  1. 10km distance to airports
  2. 300m distance to major roads
  3. no natural protection areas
  4. maximum elevation of 2000m
  5. 1000m distance to built-up areas
  6. only on suitable land cover classes

Criteria for solar implementation:

  1. only on suitable land cover classes

  2. no natural protection areas

In [6]:
# Defining Exclusion and Inclusion areas as data frames:

Exclusion_area = pd.DataFrame({"land classes" :[111, 112, 121, 122, 123, 124, 142, 511, 512, 131, 322],
                   "code value" : [1, 2, 3, 4, 5, 6, 11, 40, 41, 7, 27],}, 
                    index = ['Continuous urban fabric', 'Discontinuous urban fabric', 
                  'Industrial or commercial units', 'Road and rail networks and associated land', 
                  'Port areas', 'Airports', 'Sport and leisure facilities', 'Water courses', 'Water bodies','Mineral extraction sites','Moors and heathland'])
Exclusion_area = Exclusion_area.rename_axis(index='descriptive names')

Inclusion_area = pd.DataFrame({"land classes" :[211, 212, 213, 231, 241, 242, 243, 321, 323, 324, 333],
                   "code value" : [12, 13, 14, 18, 19, 20, 21, 26, 28, 29, 32],}, 
                    index = ['Non-irrigated arable land', 'Permanently irrigated land', 'Rice fields', 
                  'Pastures', 'Annual crops associated with permanent crops', 'Complex cultivation patterns', 
                  'Land principally occupied by agriculture with significant areas of natural vegetation', 
                  'Natural grasslands', 'Sclerophyllous vegetation', 'Transitional woodland-shrub', 
                  'Sparsely vegetated areas'])
Inclusion_area = Inclusion_area.rename_axis(index='descriptive names')

# Defining function to plot the eligible area on the country map based on excluded & included areas:
def plot_area(masked, transform, shape):
    fig, ax = plt.subplots(figsize=(5, 5))
    ax = show(masked, transform=transform, cmap='Greens', vmin=0, ax=ax)
    shape.plot(ax=ax, edgecolor='k', color='None', linewidth=2)

# Defining a function to calculate the eligible area portion of the country in percentage based on excluded and included areas:
def eligible_area(masked, excluder, shape):
    return masked.sum() * float(excluder.res**2) / shape.geometry.area.sum() * 100

####################################################################################################

# Excluding the unsuitable area for onshore wind technology:

excluder_onshore_wind = ExclusionContainer()

excluder_onshore_wind.add_geometry("ne_10m_roads.gpkg", buffer=300)                 # 300m distance to major roads
excluder_onshore_wind.add_geometry("ne_10m_airports (1).gpkg", buffer=10000)        # 10km distance to airports
excluder_onshore_wind.add_raster('WDPA_Oct2022_Public_shp-SVK.tif', crs=3035)       # no natural protection areas
excluder_onshore_wind.add_raster('PROBAV_LC100_global_v3_0_1_2019_nrt_Discrete_Classification_map.tif', codes=[1,2,3,4,5,6,11,40,41], buffer=1000, crs=3035)  #1000m distance to built-up areas

onshore_clases = [12, 13, 14, 18, 19, 20, 21, 26, 28, 29, 32]         # only on suitable land cover classes
excluder_onshore_wind.add_raster('PROBAV_LC100_global_v3_0_1_2019_nrt_Discrete_Classification_map.tif', codes=onshore_clases, crs=3035)

masked, transform = shape_availability(shape_svk, excluder_onshore_wind)

# Calculating the eligible area for onshore wind technology in percentage (%):
onshore_eligible_area = eligible_area(masked, excluder_onshore_wind, shape_svk)
print()
print(f'Almost {round(onshore_eligible_area,3)} % of Slovakia is eligible for onshore wind technology development')
print()

# Plotting the eligible area for onshore wind technology:
plot_area(masked, transform, shape_svk)
plt.text(0.5, -0.4,f'Almost {round(onshore_eligible_area,3)} % of Slovakia is eligible for onshore wind technology development', ha='center', va='center', transform=plt.gca().transAxes)

plt.savefig('on-shore wind power eligibility.png')

plt.show()
Almost 5.432 % of Slovakia is eligible for onshore wind technology development

No description has been provided for this image
In [7]:
# Excluding the unsuitable area for solar power plants:

excluder_solar = ExclusionContainer()
Solar_classes = [1, 2, 3, 4, 5, 6, 11, 18, 26, 27, 30, 31, 41]
excluder_solar.add_raster('PROBAV_LC100_global_v3_0_1_2019_nrt_Discrete_Classification_map.tif', codes=Solar_classes, crs=3035)
excluder_solar = ExclusionContainer()
excluder_solar.add_raster('WDPA_Oct2022_Public_shp-SVK.tif', crs=3035)

masked, transform = shape_availability(shape_svk, excluder_solar)

# Calculating the eligible area for solar development in percentage (%):
eligible_area_solar = eligible_area(masked, excluder_solar,shape_svk)
print()
print(f'Almost {round(eligible_area_solar,2)} % of Slovakia land is eligible for solar')
print()

# Plotting the eligible area for solar development:
plot_area(masked, transform, shape_svk)
plt.text(0.5, -0.4,f'Almost {round(eligible_area_solar,2)} % of Slovakia land is eligible for solar implementation.', ha='center', va='center', transform=plt.gca().transAxes)

plt.savefig('Solar power eligibility.png')

plt.show()
Almost 62.69 % of Slovakia land is eligible for solar

No description has been provided for this image

Secondly, another land eligibility has been performed as follows, but this time using historical weather data of ERA5 dataset & atlite library.¶

The historical weather data of Slovakia in 2017 has been downloaded into an atlite.Cutout.
A buffer of 0.25 degreesdhas also been added to the geographical bounds ofthed countr:).

In [8]:
# Define the file path for the data, being downloaded from the ERA5 dataset:
data_file_path = "era5-2017-SVK.nc"

# defining the bound limits to create a spatial extent:
minx, miny, maxx, maxy = svk_regions.total_bounds
buffer = 0.25              # A buffer of 0.25 degrees has been considered as requested.

# Creating the cutout for the year 2017:
cutout = atlite.Cutout(
    path=data_file_path,
    module="era5",
    x=slice(minx - buffer, maxx + buffer),
    y=slice(miny - buffer, maxy + buffer),
    time="2017",
)

# Calling the function cutout.prepare() & initiate the download:
cutout.prepare()

# the downloaded data is accessible here below :
cutout.data
INFO:atlite.data:Cutout already prepared.
Out[8]:
<xarray.Dataset>
Dimensions:           (x: 25, y: 10, time: 8760)
Coordinates:
  * x                 (x) float64 16.75 17.0 17.25 17.5 ... 22.25 22.5 22.75
  * y                 (y) float64 47.5 47.75 48.0 48.25 ... 49.25 49.5 49.75
  * time              (time) datetime64[ns] 2017-01-01 ... 2017-12-31T23:00:00
    lon               (x) float64 dask.array<chunksize=(25,), meta=np.ndarray>
    lat               (y) float64 dask.array<chunksize=(10,), meta=np.ndarray>
Data variables: (12/13)
    height            (y, x) float32 dask.array<chunksize=(10, 25), meta=np.ndarray>
    wnd100m           (time, y, x) float32 dask.array<chunksize=(100, 10, 25), meta=np.ndarray>
    wnd_azimuth       (time, y, x) float32 dask.array<chunksize=(100, 10, 25), meta=np.ndarray>
    roughness         (time, y, x) float32 dask.array<chunksize=(100, 10, 25), meta=np.ndarray>
    influx_toa        (time, y, x) float32 dask.array<chunksize=(100, 10, 25), meta=np.ndarray>
    influx_direct     (time, y, x) float32 dask.array<chunksize=(100, 10, 25), meta=np.ndarray>
    ...                ...
    albedo            (time, y, x) float32 dask.array<chunksize=(100, 10, 25), meta=np.ndarray>
    solar_altitude    (time, y, x) float64 dask.array<chunksize=(100, 10, 25), meta=np.ndarray>
    solar_azimuth     (time, y, x) float64 dask.array<chunksize=(100, 10, 25), meta=np.ndarray>
    temperature       (time, y, x) float32 dask.array<chunksize=(100, 10, 25), meta=np.ndarray>
    soil temperature  (time, y, x) float32 dask.array<chunksize=(100, 10, 25), meta=np.ndarray>
    runoff            (time, y, x) float32 dask.array<chunksize=(100, 10, 25), meta=np.ndarray>
Attributes:
    module:             era5
    prepared_features:  ['temperature', 'influx', 'wind', 'runoff', 'height']
    chunksize_time:     100
    Conventions:        CF-1.6
    history:            2024-01-31 15:41:18 GMT by grib_to_netcdf-2.25.1: /op...
xarray.Dataset
    • x: 25
    • y: 10
    • time: 8760
    • x
      (x)
      float64
      16.75 17.0 17.25 ... 22.5 22.75
      array([16.75, 17.  , 17.25, 17.5 , 17.75, 18.  , 18.25, 18.5 , 18.75, 19.  ,
             19.25, 19.5 , 19.75, 20.  , 20.25, 20.5 , 20.75, 21.  , 21.25, 21.5 ,
             21.75, 22.  , 22.25, 22.5 , 22.75])
    • y
      (y)
      float64
      47.5 47.75 48.0 ... 49.5 49.75
      array([47.5 , 47.75, 48.  , 48.25, 48.5 , 48.75, 49.  , 49.25, 49.5 , 49.75])
    • time
      (time)
      datetime64[ns]
      2017-01-01 ... 2017-12-31T23:00:00
      array(['2017-01-01T00:00:00.000000000', '2017-01-01T01:00:00.000000000',
             '2017-01-01T02:00:00.000000000', ..., '2017-12-31T21:00:00.000000000',
             '2017-12-31T22:00:00.000000000', '2017-12-31T23:00:00.000000000'],
            dtype='datetime64[ns]')
    • lon
      (x)
      float64
      dask.array<chunksize=(25,), meta=np.ndarray>
      Array Chunk
      Bytes 200 B 200 B
      Shape (25,) (25,)
      Dask graph 1 chunks in 2 graph layers
      Data type float64 numpy.ndarray
      25 1
    • lat
      (y)
      float64
      dask.array<chunksize=(10,), meta=np.ndarray>
      Array Chunk
      Bytes 80 B 80 B
      Shape (10,) (10,)
      Dask graph 1 chunks in 2 graph layers
      Data type float64 numpy.ndarray
      10 1
    • height
      (y, x)
      float32
      dask.array<chunksize=(10, 25), meta=np.ndarray>
      module :
      era5
      feature :
      height
      Array Chunk
      Bytes 0.98 kiB 0.98 kiB
      Shape (10, 25) (10, 25)
      Dask graph 1 chunks in 2 graph layers
      Data type float32 numpy.ndarray
      25 10
    • wnd100m
      (time, y, x)
      float32
      dask.array<chunksize=(100, 10, 25), meta=np.ndarray>
      units :
      m s**-1
      long_name :
      100 metre wind speed
      module :
      era5
      feature :
      wind
      Array Chunk
      Bytes 8.35 MiB 97.66 kiB
      Shape (8760, 10, 25) (100, 10, 25)
      Dask graph 88 chunks in 2 graph layers
      Data type float32 numpy.ndarray
      25 10 8760
    • wnd_azimuth
      (time, y, x)
      float32
      dask.array<chunksize=(100, 10, 25), meta=np.ndarray>
      units :
      m s**-1
      long_name :
      100 metre U wind component
      module :
      era5
      feature :
      wind
      Array Chunk
      Bytes 8.35 MiB 97.66 kiB
      Shape (8760, 10, 25) (100, 10, 25)
      Dask graph 88 chunks in 2 graph layers
      Data type float32 numpy.ndarray
      25 10 8760
    • roughness
      (time, y, x)
      float32
      dask.array<chunksize=(100, 10, 25), meta=np.ndarray>
      units :
      m
      long_name :
      Forecast surface roughness
      module :
      era5
      feature :
      wind
      Array Chunk
      Bytes 8.35 MiB 97.66 kiB
      Shape (8760, 10, 25) (100, 10, 25)
      Dask graph 88 chunks in 2 graph layers
      Data type float32 numpy.ndarray
      25 10 8760
    • influx_toa
      (time, y, x)
      float32
      dask.array<chunksize=(100, 10, 25), meta=np.ndarray>
      units :
      W m**-2
      module :
      era5
      feature :
      influx
      Array Chunk
      Bytes 8.35 MiB 97.66 kiB
      Shape (8760, 10, 25) (100, 10, 25)
      Dask graph 88 chunks in 2 graph layers
      Data type float32 numpy.ndarray
      25 10 8760
    • influx_direct
      (time, y, x)
      float32
      dask.array<chunksize=(100, 10, 25), meta=np.ndarray>
      units :
      W m**-2
      module :
      era5
      feature :
      influx
      Array Chunk
      Bytes 8.35 MiB 97.66 kiB
      Shape (8760, 10, 25) (100, 10, 25)
      Dask graph 88 chunks in 2 graph layers
      Data type float32 numpy.ndarray
      25 10 8760
    • influx_diffuse
      (time, y, x)
      float32
      dask.array<chunksize=(100, 10, 25), meta=np.ndarray>
      units :
      W m**-2
      module :
      era5
      feature :
      influx
      Array Chunk
      Bytes 8.35 MiB 97.66 kiB
      Shape (8760, 10, 25) (100, 10, 25)
      Dask graph 88 chunks in 2 graph layers
      Data type float32 numpy.ndarray
      25 10 8760
    • albedo
      (time, y, x)
      float32
      dask.array<chunksize=(100, 10, 25), meta=np.ndarray>
      units :
      (0 - 1)
      long_name :
      Albedo
      module :
      era5
      feature :
      influx
      Array Chunk
      Bytes 8.35 MiB 97.66 kiB
      Shape (8760, 10, 25) (100, 10, 25)
      Dask graph 88 chunks in 2 graph layers
      Data type float32 numpy.ndarray
      25 10 8760
    • solar_altitude
      (time, y, x)
      float64
      dask.array<chunksize=(100, 10, 25), meta=np.ndarray>
      time shift :
      -1 days +23:30:00
      units :
      rad
      module :
      era5
      feature :
      influx
      Array Chunk
      Bytes 16.71 MiB 195.31 kiB
      Shape (8760, 10, 25) (100, 10, 25)
      Dask graph 88 chunks in 2 graph layers
      Data type float64 numpy.ndarray
      25 10 8760
    • solar_azimuth
      (time, y, x)
      float64
      dask.array<chunksize=(100, 10, 25), meta=np.ndarray>
      time shift :
      -1 days +23:30:00
      units :
      rad
      module :
      era5
      feature :
      influx
      Array Chunk
      Bytes 16.71 MiB 195.31 kiB
      Shape (8760, 10, 25) (100, 10, 25)
      Dask graph 88 chunks in 2 graph layers
      Data type float64 numpy.ndarray
      25 10 8760
    • temperature
      (time, y, x)
      float32
      dask.array<chunksize=(100, 10, 25), meta=np.ndarray>
      units :
      K
      long_name :
      2 metre temperature
      module :
      era5
      feature :
      temperature
      Array Chunk
      Bytes 8.35 MiB 97.66 kiB
      Shape (8760, 10, 25) (100, 10, 25)
      Dask graph 88 chunks in 2 graph layers
      Data type float32 numpy.ndarray
      25 10 8760
    • soil temperature
      (time, y, x)
      float32
      dask.array<chunksize=(100, 10, 25), meta=np.ndarray>
      units :
      K
      long_name :
      Soil temperature level 4
      module :
      era5
      feature :
      temperature
      Array Chunk
      Bytes 8.35 MiB 97.66 kiB
      Shape (8760, 10, 25) (100, 10, 25)
      Dask graph 88 chunks in 2 graph layers
      Data type float32 numpy.ndarray
      25 10 8760
    • runoff
      (time, y, x)
      float32
      dask.array<chunksize=(100, 10, 25), meta=np.ndarray>
      units :
      m
      long_name :
      Runoff
      module :
      era5
      feature :
      runoff
      Array Chunk
      Bytes 8.35 MiB 97.66 kiB
      Shape (8760, 10, 25) (100, 10, 25)
      Dask graph 88 chunks in 2 graph layers
      Data type float32 numpy.ndarray
      25 10 8760
    • x
      PandasIndex
      PandasIndex(Index([16.75,  17.0, 17.25,  17.5, 17.75,  18.0, 18.25,  18.5, 18.75,  19.0,
             19.25,  19.5, 19.75,  20.0, 20.25,  20.5, 20.75,  21.0, 21.25,  21.5,
             21.75,  22.0, 22.25,  22.5, 22.75],
            dtype='float64', name='x'))
    • y
      PandasIndex
      PandasIndex(Index([47.5, 47.75, 48.0, 48.25, 48.5, 48.75, 49.0, 49.25, 49.5, 49.75], dtype='float64', name='y'))
    • time
      PandasIndex
      PandasIndex(DatetimeIndex(['2017-01-01 00:00:00', '2017-01-01 01:00:00',
                     '2017-01-01 02:00:00', '2017-01-01 03:00:00',
                     '2017-01-01 04:00:00', '2017-01-01 05:00:00',
                     '2017-01-01 06:00:00', '2017-01-01 07:00:00',
                     '2017-01-01 08:00:00', '2017-01-01 09:00:00',
                     ...
                     '2017-12-31 14:00:00', '2017-12-31 15:00:00',
                     '2017-12-31 16:00:00', '2017-12-31 17:00:00',
                     '2017-12-31 18:00:00', '2017-12-31 19:00:00',
                     '2017-12-31 20:00:00', '2017-12-31 21:00:00',
                     '2017-12-31 22:00:00', '2017-12-31 23:00:00'],
                    dtype='datetime64[ns]', name='time', length=8760, freq=None))
  • module :
    era5
    prepared_features :
    ['temperature', 'influx', 'wind', 'runoff', 'height']
    chunksize_time :
    100
    Conventions :
    CF-1.6
    history :
    2024-01-31 15:41:18 GMT by grib_to_netcdf-2.25.1: /opt/ecmwf/mars-client/bin/grib_to_netcdf.bin -S param -o /cache/data5/adaptor.mars.internal-1706715678.4543219-7391-1-122bc265-b6ee-4f49-8125-b6c519ca7a5e.nc /cache/tmp/122bc265-b6ee-4f49-8125-b6c519ca7a5e-adaptor.mars.internal-1706715663.8918104-7391-1-tmp.grib
In [9]:
# Which data is contained in the cutout:
cutout
Out[9]:
<Cutout "era5-2017-SVK">
 x = 16.75 ⟷ 22.75, dx = 0.25
 y = 47.50 ⟷ 49.75, dy = 0.25
 time = 2017-01-01 ⟷ 2017-12-31, dt = H
 module = era5
 prepared_features = ['height', 'wind', 'influx', 'temperature', 'runoff']
In [10]:
# Included weather variables list:

cutout.prepared_features
Out[10]:
module  feature    
era5    height                   height
        wind                    wnd100m
        wind                wnd_azimuth
        wind                  roughness
        influx               influx_toa
        influx            influx_direct
        influx           influx_diffuse
        influx                   albedo
        influx           solar_altitude
        influx            solar_azimuth
        temperature         temperature
        temperature    soil temperature
        runoff                   runoff
dtype: object
In [11]:
# Creating availability matrixes for solar:
?cutout.availabilitymatrix
solar_AM = cutout.availabilitymatrix(shape_svk, excluder_solar)
#display(solar_AM)

# Creating availability matrixes for onshore wind:
onshore_AM = cutout.availabilitymatrix(shape_svk, excluder_onshore_wind)
#display(onshore_AM)

cap_per_sqkm = 3          # [MW/km2]  For both wind and solar, a deployment density of 3 MW/km2 has been assumed.

# Calculating the area :
country_area = cutout.grid.set_index(["y", "x"]).to_crs(3035).area / 1e6
country_area = xr.DataArray(country_area, dims=("spatial"))
#display(country_area)
Compute availability matrix: 100%|███████████████████████████████████████████████| 8/8 [00:01<00:00,  4.56 gridcells/s]
Compute availability matrix: 100%|███████████████████████████████████████████████| 8/8 [01:43<00:00, 12.97s/ gridcells]
Signature:
cutout.availabilitymatrix(
    shapes,
    excluder,
    nprocesses=None,
    disable_progressbar=False,
)
Docstring:
Compute the eligible share within cutout cells in the overlap with shapes.

For parallel calculation (nprocesses not None) the excluder must not be
initialized and all raster references must be strings. Otherwise processes
are colliding when reading from one common rasterio.DatasetReader.

Parameters
----------
cutout : atlite.Cutout
    Cutout which the availability matrix is aligned to.
shapes : geopandas.Series/geopandas.DataFrame
    Geometries for which the availabilities are calculated.
excluder : atlite.gis.ExclusionContainer
    Container of all meta data or objects which to exclude, i.e.
    rasters and geometries.
nprocesses : int, optional
    Number of processes to use for calculating the matrix. The paralle-
    lization can heavily boost the calculation speed. The default is None.
disable_progressbar: bool, optional
    Disable the progressbar if nprocesses is not None. Then the `map`
    function instead of the `imap` function is used for the multiprocessing
    pool. This speeds up the calculation.

Returns
-------
availabilities : xr.DataArray
    DataArray of shape (|shapes|, |y|, |x|) containing all the eligible
    share of cutout cell (x,y) in the overlap with shape i.

Notes
-----
The rasterio (or GDAL) average downsampling returns different results
dependent on how the target raster (the cutout raster) is spanned.
Either it is spanned from the top left going downwards,
e.g. Affine(0.25, 0, 0, 0, -0.25, 50), or starting in the
lower left corner and going up, e.g. Affine(0.25, 0, 0, 0, 0.25, 50).
Here we stick to the top down version which is why we use
`cutout.transform_r` and flipping the y-axis in the end.
File:      c:\users\mahtab\anaconda3\envs\esm-2023\lib\site-packages\atlite\gis.py
Type:      method
In [12]:
# Creating capacity matrix for onshore wind:
capacity_matrix_onshore = onshore_AM.stack(spatial=["y", "x"]) * country_area * cap_per_sqkm
#display(capacity_matrix_onshore)

# Wind profile:
wind_onshore = cutout.wind(
    atlite.windturbines.Vestas_V112_3MW,
    matrix=capacity_matrix_onshore,
    index= svk_regions.index,
    per_unit=True,
)

svk_wind_onshore_2017_potential_time_series = wind_onshore.to_pandas()
display(svk_wind_onshore_2017_potential_time_series.head(15))
INFO:atlite.convert:Convert and aggregate 'wind'.
NAME_1 Banskobystrický Bratislavský Košický Nitriansky Prešovský Trenčiansky Trnavský Žilinský
time
2017-01-01 00:00:00 0.000003 0.000060 0.007087 0.000000e+00 0.030756 0.000335 0.000001 0.013306
2017-01-01 01:00:00 0.000000 0.000005 0.007612 0.000000e+00 0.034854 0.000475 0.000040 0.014443
2017-01-01 02:00:00 0.000000 0.000066 0.007142 1.130307e-06 0.039974 0.000691 0.000331 0.015357
2017-01-01 03:00:00 0.000000 0.000171 0.003622 5.924064e-07 0.041654 0.001257 0.000336 0.017744
2017-01-01 04:00:00 0.000000 0.000144 0.003500 0.000000e+00 0.047368 0.001554 0.000341 0.020749
2017-01-01 05:00:00 0.000000 0.000156 0.003398 0.000000e+00 0.054182 0.002014 0.000318 0.027655
2017-01-01 06:00:00 0.000000 0.000167 0.003533 4.355234e-08 0.062805 0.002464 0.000342 0.035510
2017-01-01 07:00:00 0.000000 0.000406 0.004653 6.033786e-07 0.067074 0.002595 0.000470 0.042285
2017-01-01 08:00:00 0.000000 0.000687 0.004192 2.199573e-06 0.066282 0.002300 0.001231 0.044154
2017-01-01 09:00:00 0.000000 0.000343 0.005305 1.851060e-07 0.058913 0.001037 0.000373 0.044615
2017-01-01 10:00:00 0.000019 0.000218 0.005946 1.088100e-06 0.053568 0.000092 0.000304 0.027705
2017-01-01 11:00:00 0.000000 0.000181 0.007191 1.481370e-06 0.049395 0.000112 0.000274 0.025161
2017-01-01 12:00:00 0.000000 0.000299 0.008346 1.112223e-04 0.067957 0.000216 0.000336 0.031990
2017-01-01 13:00:00 0.000000 0.000523 0.008210 1.796761e-04 0.085420 0.000869 0.000486 0.044031
2017-01-01 14:00:00 0.000106 0.001866 0.008984 3.887791e-03 0.107708 0.005368 0.007772 0.060552
In [13]:
# Plotting the on-shore wind potential time series for Jan-2017: 

plt.figure(figsize=(15, 6))  

region_colors = {
    'Banskobystrický': 'b',  
    'Bratislavský': 'g',  
    'Košický': 'r',  
    'Nitriansky': 'y',  
    'Prešovský': 'brown',  
    'Trenčiansky': 'pink',  
    'Trnavský': 'k', 
    'Žilinský': 'c', 
}

for region, color in region_colors.items():
    svk_wind_onshore_2017_potential_time_series.loc['2017-01', region].plot(color=color, label=region)

plt.ylabel('On-shore Wind Potential [p.u.]')  
plt.ylim(0, 1)
plt.title('On-shore Wind Potential [p.u.] of all regions in the January, 2017')
plt.savefig('On-shore Wind Potential-January.png')
plt.legend()
plt.grid(True)
plt.show()


################################################

# Plotting the on-shore wind potential time series: ( Banskobystrický)

plt.figure(figsize=(12, 6))  
svk_wind_onshore_2017_potential_time_series.iloc[:,0].plot(color='red')
plt.ylabel('On-shore Wind Potential [p.u.]')  
plt.ylim(0, 1)
plt.title('On-shore Wind Potential [p.u.] for Banskobystrický region in the country of Slovokia')
plt.savefig('On-shore Wind Potential-Banskobystrický.png')

plt.show()

##################################################
# Plotting the on-shore wind potential time series: ( Bratislavský)

plt.figure(figsize=(12, 6))  
svk_wind_onshore_2017_potential_time_series.iloc[:,1].plot(color='red')
plt.ylabel('On-shore Wind Potential [p.u.]')  
plt.ylim(0, 1)
plt.title('On-shore Wind Potential [p.u.] for Bratislavský region in the country of Slovokia')
plt.savefig('On-shore Wind Potential-Bratislavský.png')

plt.show()

################################################

# Plotting the on-shore wind potential time series: ( Košický)

plt.figure(figsize=(12, 6))  
svk_wind_onshore_2017_potential_time_series.iloc[:,2].plot(color='green')
plt.ylabel('On-shore Wind Potential [p.u.]')  
plt.ylim(0, 1)
plt.title('On-shore Wind Potential [p.u.] for Košický region in the country of Slovokia')
plt.savefig('On-shore Wind Potential-Košický.png')

plt.show()


################################################

# Plotting the on-shore wind potential time series: ( Nitriansky)

plt.figure(figsize=(12, 6))  
svk_wind_onshore_2017_potential_time_series.iloc[:,3].plot(color='grey')
plt.ylabel('On-shore Wind Potential [p.u.]')  
plt.ylim(0, 1)
plt.title('On-shore Wind Potential [p.u.] for Nitriansky region in the country of Slovokia')
plt.savefig('On-shore Wind Potential-Nitriansky.png')

plt.show()
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
In [14]:
# Creating capacity matrix for solar:
capacity_matrix_solar = solar_AM.stack(spatial=['y', 'x']) * country_area * cap_per_sqkm
#display(capacity_matrix_solar)

# Solar PV profile:
pv = cutout.pv(
    panel=atlite.solarpanels.CdTe,
    matrix=capacity_matrix_solar,
    orientation="latitude_optimal",
    index= svk_regions.index,
    per_unit=True,
)

svk_pv_2017_potential_time_series = pv.to_pandas()
svk_pv_2017_potential_time_series.tail(15)
INFO:atlite.convert:Convert and aggregate 'pv'.
[########################################] | 100% Completed | 5.36 s
Out[14]:
NAME_1 Banskobystrický Bratislavský Košický Nitriansky Prešovský Trenčiansky Trnavský Žilinský
time
2017-12-31 09:00:00 0.028286 0.013923 0.074632 0.020029 0.142586 0.017215 0.014964 0.035910
2017-12-31 10:00:00 0.049915 0.022226 0.118199 0.042010 0.244426 0.039485 0.027253 0.065328
2017-12-31 11:00:00 0.156638 0.026005 0.158670 0.046129 0.304047 0.045949 0.020684 0.143634
2017-12-31 12:00:00 0.085264 0.050059 0.127709 0.042245 0.222590 0.027803 0.034511 0.067804
2017-12-31 13:00:00 0.050969 0.161593 0.068615 0.061506 0.163986 0.024984 0.110701 0.047229
2017-12-31 14:00:00 0.038435 0.237795 0.022509 0.117132 0.056720 0.030638 0.184676 0.016817
2017-12-31 15:00:00 0.012403 0.129751 0.005537 0.069491 0.010469 0.028245 0.108905 0.003444
2017-12-31 16:00:00 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
2017-12-31 17:00:00 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
2017-12-31 18:00:00 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
2017-12-31 19:00:00 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
2017-12-31 20:00:00 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
2017-12-31 21:00:00 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
2017-12-31 22:00:00 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
2017-12-31 23:00:00 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
In [15]:
# Plotting the solar potential time series for Jan-2017: 
plt.figure(figsize=(15, 6))  

for region, color in region_colors.items():
    pv.loc['2017-01', region].plot(color=color, label=region)

plt.ylabel('On-shore Wind Potential [p.u.]')  
plt.ylim(0, 1)
plt.title('Solar Potential [p.u.] of all regions in the January, 2017')
plt.savefig('Solar Potential-January.png')
plt.legend()
plt.grid(True)
plt.show()

##########################################################


# Plotting the solar potential time series: (only Banskobystrický)

plt.figure(figsize=(12, 6)) 
pv.to_pandas().iloc[:,0].plot(color='blue')
plt.ylabel('Solar Potential [p.u.]')  
plt.ylim(0, 1)

plt.title('Solar Potential [p.u.] for Banskobystrický region in the country of Slovokia')
plt.savefig('Solar Potential-Banskobystrický.png')

plt.show()

##################################################

# Plotting the solar potential time series: (only Bratislavský)

plt.figure(figsize=(12, 6)) 
pv.to_pandas().iloc[:,1].plot(color='red')
plt.ylabel('Solar Potential [p.u.]')  
plt.ylim(0, 1)

plt.title('Solar Potential [p.u.] for Bratislavský region in the country of Slovokia')
plt.savefig('Solar Potential-Bratislavský.png')

plt.show()


######################################################


# Plotting the solar potential time series: (only Košický)

plt.figure(figsize=(12, 6)) 
pv.to_pandas().iloc[:,2].plot(color='green')
plt.ylabel('Solar Potential [p.u.]')  
plt.ylim(0, 1)

plt.title('Solar Potential [p.u.] for Košický region in the country of Slovokia')
plt.savefig('Solar Potential-Košický.png')

plt.show()

################################################################

# Plotting the solar potential time series: (only Nitriansky)

plt.figure(figsize=(12, 6)) 
pv.to_pandas().iloc[:,3].plot(color='g')
plt.ylabel('Solar Potential [p.u.]')  
plt.ylim(0, 1)

plt.title('Solar Potential [p.u.] for Nitriansky region in the country of Slovokia')
plt.savefig('Solar Potential-Nitriansky.png')

plt.show()
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image

Building the Model:¶

Below, a PyPSA model is being created that minimizes total annual system costs with the specified characteristics, according to the assignment explanations:¶

In [48]:
from matplotlib.cm import ScalarMappable

## Building the Pypsa Model to minimize total annual system costs:

#display(gdf_svk.head(3))


n = pypsa.Network()       # defining a network

start = '2017-01-01 00:00:00'
end = '2017-12-31 23:00:00'

timestamps = pd.date_range(start=start, end=end, freq='H')

n.set_snapshots(timestamps)


# Defining the Buses:

regions_points = svk_regions.representative_point()

for index, row in svk_regions.iterrows():                    # A bus (Node) is defined for each region based on its geometry (latitude & Longitude)
    bus_name = index
    if bus_name in n.buses.index:
        n.remove("Bus", bus_name)
    region_point = regions_points.loc[index]  
    n.add("Bus", bus_name, x=region_point.x, y=region_point.y)

display(n.buses)

# !!!! FOR PRESENTATION: !!!
fig, ax = plt.subplots(figsize=(13, 9))
svk_regions.plot(ax=ax, facecolor='lightblue', edgecolor='black', linewidth=1.5)
regions_points.plot(ax=ax, color="black", markersize=20);
for bus_name, bus in n.buses.iterrows():
    ax.annotate(bus_name, (bus['x'], bus['y']), textcoords="offset points", xytext=(3,10), ha='center')
plt.title('regions of Slovokia as Nodes in the Model')
plt.savefig('Nodes.png')

plt.show()

##################################################################
# !!!! FOR PRESENTATION: !!!
fig, ax = plt.subplots(figsize=(13, 9))
svk_regions.plot(ax=ax, facecolor='lightblue', edgecolor='black', linewidth=1.5)
regions_points.plot(ax=ax, color="lightblue", markersize=0.5);


for bus_name, bus in n.buses.iterrows():
    ax.annotate(bus_name, (bus['x'], bus['y']), textcoords="offset points", xytext=(3,10), ha='center')

power_plants_scatter = gdf_svk.plot(
    ax=ax ,
    column = gdf_svk['primary_fuel'],
    markersize = gdf_svk['capacity_mw'] / 1e1 *4,
    legend= True,
    cmap='viridis',
)

plt.title('Power plants Capacities per region in Slovokia')
plt.savefig('Power plants capacities.png')

plt.show()
attribute v_nom type x y carrier unit v_mag_pu_set v_mag_pu_min v_mag_pu_max control sub_network
Bus
Banskobystrický 1.0 19.424031 48.499321 AC 1.0 0.0 inf PQ
Bratislavský 1.0 17.175805 48.331983 AC 1.0 0.0 inf PQ
Košický 1.0 21.292701 48.676239 AC 1.0 0.0 inf PQ
Nitriansky 1.0 18.336016 48.223509 AC 1.0 0.0 inf PQ
Prešovský 1.0 21.219622 49.114946 AC 1.0 0.0 inf PQ
Trenčiansky 1.0 18.215575 48.902964 AC 1.0 0.0 inf PQ
Trnavský 1.0 17.645627 48.323494 AC 1.0 0.0 inf PQ
Žilinský 1.0 19.196299 49.177544 AC 1.0 0.0 inf PQ
No description has been provided for this image
No description has been provided for this image
In [17]:
year = 2040    # A year for cost projections has been selected
url = f"https://raw.githubusercontent.com/PyPSA/technology-data/master/outputs/costs_{year}.csv"
tech_properties = pd.read_csv(url)
#display(tech_properties.head(10))

tech_properties_svk = tech_properties[
    (tech_properties['technology'] == 'nuclear') |
    (tech_properties['technology'] == 'CCGT') |
    (tech_properties['technology'] == 'coal') |
    (tech_properties['technology'] == 'hydro') |
    (tech_properties['technology'] == 'solar') |
    (tech_properties['technology'] == 'gas') |
    (tech_properties['technology'] == 'electrolysis') |
    (tech_properties['technology'] == 'fuel cell') |
    (tech_properties['technology'] == 'hydrogen storage underground') |
    (tech_properties['technology'] == 'battery inverter') |
    (tech_properties['technology'] == 'battery storage')
    
]
tech_properties_svk['technology'] = tech_properties_svk['technology'].replace('gas', 'CCGT')
display(tech_properties_svk)
technology parameter value unit source further description
29 CCGT FOM 3.3006 %/year Danish Energy Agency, technology_data_for_el_a... 05 Gas turb. CC, steam extract.: Fixed O&M
30 CCGT VOM 4.1000 EUR/MWh Danish Energy Agency, technology_data_for_el_a... 05 Gas turb. CC, steam extract.: Variable O&M
31 CCGT c_b 2.1000 50oC/100oC Danish Energy Agency, technology_data_for_el_a... 05 Gas turb. CC, steam extract.: Cb coefficient
32 CCGT c_v 0.1500 50oC/100oC Danish Energy Agency, technology_data_for_el_a... 05 Gas turb. CC, steam extract.: Cv coefficient
33 CCGT efficiency 0.5900 per unit Danish Energy Agency, technology_data_for_el_a... 05 Gas turb. CC, steam extract.: Electricity ...
34 CCGT investment 815.0000 EUR/kW Danish Energy Agency, technology_data_for_el_a... 05 Gas turb. CC, steam extract.: Nominal inve...
35 CCGT lifetime 25.0000 years Danish Energy Agency, technology_data_for_el_a... 05 Gas turb. CC, steam extract.: Technical li...
412 battery inverter FOM 0.5400 %/year Danish Energy Agency, technology_data_catalogu... : Fixed O&M
413 battery inverter efficiency 0.9600 per unit Danish Energy Agency, technology_data_catalogu... : Round trip efficiency DC
414 battery inverter investment 100.0000 EUR/kW Danish Energy Agency, technology_data_catalogu... : Output capacity expansion cost investment
415 battery inverter lifetime 10.0000 years Danish Energy Agency, technology_data_catalogu... : Technical lifetime
416 battery storage investment 94.0000 EUR/kWh Danish Energy Agency, technology_data_catalogu... : Energy storage expansion cost investment
417 battery storage lifetime 30.0000 years Danish Energy Agency, technology_data_catalogu... : Technical lifetime
579 coal CO2 intensity 0.3361 tCO2/MWh_th Entwicklung der spezifischen Kohlendioxid-Emis... NaN
580 coal FOM 1.3100 %/year Lazard's levelized cost of energy analysis - v... Calculated based on average of listed range, i...
581 coal VOM 3.3278 EUR/MWh_e Lazard's levelized cost of energy analysis - v... Calculated based on average of listed range, i...
582 coal efficiency 0.3300 p.u. Lazard's levelized cost of energy analysis - v... Calculated based on average of listed range, i...
583 coal fuel 9.2743 EUR/MWh_th DIW (2013): Current and propsective costs of e... Based on IEA 2011 data, 99 USD/t.
584 coal investment 3905.3074 EUR/kW_e Lazard's levelized cost of energy analysis - v... Higher costs include coal plants with CCS, but...
585 coal lifetime 40.0000 years Lazard's levelized cost of energy analysis - v... NaN
707 electrolysis FOM 4.0000 %/year Danish Energy Agency, data_sheets_for_renewabl... 86 AEC 100 MW: Fixed O&M
708 electrolysis efficiency 0.6532 per unit Danish Energy Agency, data_sheets_for_renewabl... 86 AEC 100 MW: Hydrogen Output
709 electrolysis efficiency-heat 0.1849 per unit Danish Energy Agency, data_sheets_for_renewabl... 86 AEC 100 MW: - hereof recoverable for dist...
710 electrolysis investment 384.9356 EUR/kW_e Danish Energy Agency, data_sheets_for_renewabl... 86 AEC 100 MW: Specific investment
711 electrolysis lifetime 25.0000 years Danish Energy Agency, data_sheets_for_renewabl... 86 AEC 100 MW: Technical lifetime
712 fuel cell FOM 5.0000 %/year Danish Energy Agency, technology_data_for_el_a... 12 LT-PEMFC CHP: Fixed O&M
713 fuel cell c_b 1.2500 50oC/100oC Danish Energy Agency, technology_data_for_el_a... 12 LT-PEMFC CHP: Cb coefficient
714 fuel cell efficiency 0.5000 per unit Danish Energy Agency, technology_data_for_el_a... 12 LT-PEMFC CHP: Electricity efficiency, annu...
715 fuel cell investment 950.0000 EUR/kW_e Danish Energy Agency, technology_data_for_el_a... 12 LT-PEMFC CHP: Nominal investment
716 fuel cell lifetime 10.0000 years Danish Energy Agency, technology_data_for_el_a... 12 LT-PEMFC CHP: Technical lifetime
717 CCGT CO2 intensity 0.1980 tCO2/MWh_th Stoichiometric calculation with 50 GJ/t CH4 NaN
718 CCGT fuel 23.8481 EUR/MWh_th DIW (2013): Current and propsective costs of e... Based on IEA 2011 data.
745 hydro FOM 1.0000 %/year DIW DataDoc http://hdl.handle.net/10419/80348 from old pypsa cost assumptions
746 hydro efficiency 0.9000 per unit DIW DataDoc http://hdl.handle.net/10419/80348 from old pypsa cost assumptions
747 hydro investment 2208.1616 EUR/kWel DIW DataDoc http://hdl.handle.net/10419/80348 from old pypsa cost assumptions
748 hydro lifetime 80.0000 years IEA2010 from old pypsa cost assumptions
760 hydrogen storage underground FOM 0.0000 %/year Danish Energy Agency, technology_data_catalogu... 151c Hydrogen Storage - Caverns: Fixed O&M
761 hydrogen storage underground VOM 0.0000 EUR/MWh Danish Energy Agency, technology_data_catalogu... 151c Hydrogen Storage - Caverns: Variable O&M
762 hydrogen storage underground investment 1.5000 EUR/kWh Danish Energy Agency, technology_data_catalogu... 151c Hydrogen Storage - Caverns: Specific inv...
763 hydrogen storage underground lifetime 100.0000 years Danish Energy Agency, technology_data_catalogu... 151c Hydrogen Storage - Caverns: Technical li...
814 nuclear FOM 1.2700 %/year Lazard's levelized cost of energy analysis - v... U.S. specific costs including newly commission...
815 nuclear VOM 3.6188 EUR/MWh_e Lazard's levelized cost of energy analysis - v... U.S. specific costs including newly commission...
816 nuclear efficiency 0.3260 p.u. Lazard's levelized cost of energy analysis - v... Based on heat rate of 10.45 MMBtu/MWh_e and 3....
817 nuclear fuel 3.3122 EUR/MWh_th DIW (2013): Current and propsective costs of e... Based on IEA 2011 data.
818 nuclear investment 8769.6136 EUR/kW_e Lazard's levelized cost of energy analysis - v... U.S. specific costs including newly commission...
819 nuclear lifetime 40.0000 years Lazard's levelized cost of energy analysis - v... NaN
858 solar FOM 2.0400 %/year Calculated. See 'further description'. Mixed investment costs based on average of 50%...
859 solar VOM 0.0100 EUR/MWhel RES costs made up to fix curtailment order from old pypsa cost assumptions
860 solar investment 407.8706 EUR/kW_e Calculated. See 'further description'. Mixed investment costs based on average of 50%...
861 solar lifetime 40.0000 years Calculated. See 'further description'. Mixed investment costs based on average of 50%...
In [18]:
# Dataframe of the technologies:

df_CCGT = tech_properties_svk[tech_properties_svk.iloc[:, 0] == 'CCGT']
df_CCGT.set_index('parameter', inplace=True)
#display(df_CCGT)

df_solar = tech_properties_svk[tech_properties_svk.iloc[:, 0] == 'solar']
df_solar.set_index('parameter', inplace=True)
df_solar.loc['fuel'] = ['solar', 0, None, None, None]
df_solar.loc['efficiency'] = ['solar', 1, None, None, None]
#display(df_solar)

df_nuclear = tech_properties_svk[tech_properties_svk.iloc[:, 0] == 'nuclear']
df_nuclear.set_index('parameter', inplace=True)
#display(df_nuclear)

df_hydro = tech_properties_svk[tech_properties_svk.iloc[:, 0] == 'hydro']
df_hydro.set_index('parameter', inplace=True)
df_hydro.loc['fuel'] = ['hydro', 0, None, None, None]
df_hydro.loc['VOM'] = ['hydro', 0, None, None, None]
#display(df_hydro)

df_coal = tech_properties_svk[tech_properties_svk.iloc[:, 0] == 'coal']
df_coal.set_index('parameter', inplace=True)
#display(df_coal)

df_electrolysis = tech_properties_svk[tech_properties_svk.iloc[:, 0] == 'electrolysis']
df_electrolysis.set_index('parameter', inplace=True)

df_fuel_cell = tech_properties_svk[tech_properties_svk.iloc[:, 0] == 'fuel cell']
df_fuel_cell.set_index('parameter', inplace=True)

df_hydrogen_storage_underground = tech_properties_svk[tech_properties_svk.iloc[:, 0] == 'hydrogen storage underground']
df_hydrogen_storage_underground.set_index('parameter', inplace=True)
display(df_electrolysis)

df_battery_inverter = tech_properties_svk[tech_properties_svk.iloc[:, 0] == 'battery inverter']
df_battery_inverter.set_index('parameter', inplace=True)

df_battery_storage = tech_properties_svk[tech_properties_svk.iloc[:, 0] == 'battery storage']
df_battery_storage.set_index('parameter', inplace=True)
technology value unit source further description
parameter
FOM electrolysis 4.0000 %/year Danish Energy Agency, data_sheets_for_renewabl... 86 AEC 100 MW: Fixed O&M
efficiency electrolysis 0.6532 per unit Danish Energy Agency, data_sheets_for_renewabl... 86 AEC 100 MW: Hydrogen Output
efficiency-heat electrolysis 0.1849 per unit Danish Energy Agency, data_sheets_for_renewabl... 86 AEC 100 MW: - hereof recoverable for dist...
investment electrolysis 384.9356 EUR/kW_e Danish Energy Agency, data_sheets_for_renewabl... 86 AEC 100 MW: Specific investment
lifetime electrolysis 25.0000 years Danish Energy Agency, data_sheets_for_renewabl... 86 AEC 100 MW: Technical lifetime
In [19]:
# Calculating the marginal costs of the technologies:

marginalcost_Coal = df_coal.loc['fuel', 'value'] / df_coal.loc['efficiency', 'value'] + df_coal.loc['VOM', 'value']
marginalcost_Solar = df_solar.loc['fuel', 'value'] / df_solar.loc['efficiency', 'value'] + df_solar.loc['VOM', 'value']
marginalcost_Nuclear = df_nuclear.loc['fuel', 'value'] / df_nuclear.loc['efficiency', 'value'] + df_nuclear.loc['VOM', 'value']
marginalcost_Gas = df_CCGT.loc['fuel', 'value'] / df_CCGT.loc['efficiency', 'value'] + df_CCGT.loc['VOM', 'value']
marginalcost_Hydro = df_hydro.loc['fuel', 'value'] / df_hydro.loc['efficiency', 'value'] + df_hydro.loc['VOM', 'value']
marginalcost_Wind = 0


marginal_costs ={ 'Coal' : marginalcost_Coal,
                 'Solar' : marginalcost_Solar,
                 'Nuclear': marginalcost_Nuclear,
                 'Gas': marginalcost_Gas,
                 'Hydro' : marginalcost_Hydro,
                 'Wind' : marginalcost_Wind,
                }
In [20]:
# Calculating the capital costs of the technologies:

r = 0.07

annuity_coal = (r / (1 - (1 + r) ** -df_coal.loc['lifetime', 'value'])) * (df_coal.loc['investment', 'value']) * 1000
annuity_solar = (r / (1 - (1 + r) ** -df_solar.loc['lifetime', 'value'])) * (df_solar.loc['investment', 'value']) * 1000
annuity_hydro = (r / (1 - (1 + r) ** -df_hydro.loc['lifetime', 'value'])) * (df_hydro.loc['investment', 'value']) * 1000
annuity_nuclear = (r / (1 - (1 + r) ** -df_nuclear.loc['lifetime', 'value'])) * (df_nuclear.loc['investment', 'value']) * 1000
annuity_CCGT = (r / (1 - (1 + r) ** -df_CCGT.loc['lifetime', 'value'])) * (df_CCGT.loc['investment', 'value']) * 1000
annuity_electrolysis = (r / (1 - (1 + r) ** -df_electrolysis.loc['lifetime', 'value'])) * (df_electrolysis.loc['investment', 'value']) * 1000
annuity_fuel_cell = (r / (1 - (1 + r) ** -df_fuel_cell.loc['lifetime', 'value'])) * (df_fuel_cell.loc['investment', 'value']) * 1000
annuity_hydrogen_storage_underground = (r / (1 - (1 + r) ** -df_hydrogen_storage_underground.loc['lifetime', 'value'])) * (df_hydrogen_storage_underground.loc['investment', 'value']) * 1000
annuity_battery_inverter = (r / (1 - (1 + r) ** -df_battery_inverter.loc['lifetime', 'value'])) * (df_battery_inverter.loc['investment', 'value']) * 1000
annuity_battery_storage = (r / (1 - (1 + r) ** -df_battery_storage.loc['lifetime', 'value'])) * (df_battery_storage.loc['investment', 'value']) * 1000

FOM_coal = (df_coal.loc['FOM', 'value'] / 100) * df_coal.loc['investment', 'value'] *1000
FOM_solar = (df_solar.loc['FOM', 'value'] / 100) * df_solar.loc['investment', 'value'] *1000
FOM_nuclear = (df_nuclear.loc['FOM', 'value'] / 100) * df_nuclear.loc['investment', 'value'] *1000
FOM_hydro = (df_hydro.loc['FOM', 'value'] / 100) * df_hydro.loc['investment', 'value'] *1000
FOM_CCGT = (df_CCGT.loc['FOM', 'value'] / 100) * df_CCGT.loc['investment', 'value'] *1000
FOM_electrolysis = (df_electrolysis.loc['FOM', 'value'] / 100) * df_electrolysis.loc['investment', 'value'] *1000
FOM_fuel_cell = (df_fuel_cell.loc['FOM', 'value'] / 100) * df_fuel_cell.loc['investment', 'value'] *1000
FOM_battery_inverter = (df_battery_inverter.loc['FOM', 'value'] / 100) * df_battery_inverter.loc['investment', 'value'] *1000


capitalcost_coal = annuity_coal + FOM_coal
capitalcost_solar = annuity_solar + FOM_solar
capitalcost_nuclear = annuity_nuclear + FOM_nuclear
capitalcost_hydro = annuity_hydro + FOM_hydro
capitalcost_CCGT = annuity_CCGT + FOM_CCGT
capitalcost_electrolysis= annuity_electrolysis + FOM_electrolysis
capitalcost_fuel_cell= annuity_fuel_cell + FOM_fuel_cell
capitalcost_hydrogen_storage_underground = annuity_hydrogen_storage_underground 
capitalcost_battery_storage = annuity_battery_storage
capitalcost_battery_inverter= annuity_battery_inverter + FOM_battery_inverter

display(capitalcost_battery_storage)
7575.121930044452
In [21]:
## Defining the Generators:

power_plants_df_with_regions = gpd.sjoin(gdf_svk, svk_regions, how="left", op="within")  # spatial join to specify in which region each power plant is located
#display(power_plants_df_with_regions.head(5))

aggregated_power_plants = power_plants_df_with_regions.groupby(['primary_fuel', 'index_right']).first()  # grouping the power plant data based on technology and region

#display(aggregated_power_plants)

# Adding the Carriers & Corresponding emissions:
emissions = {
    'Coal': 0.34,
    'Gas': 0.2,
    'Hydro': 0,
    'Nuclear': 0,
    'Solar': 0,
    'Wind': 0
}

n.madd(
    "Carrier",
    ["Coal", "Gas", "Hydro", "Wind", "Solar", "Nuclear"],
    co2_emissions=emissions,
)


# Adding the generators to the network

for index, row in aggregated_power_plants.iterrows():
    generator_name = f"{index[0]}_{index[1]}"       # Specifying a name for each generator
    
    if generator_name in n.generators.index:        # Removing already created generators to avoid errors regarding duplication of generators name
        n.remove("Generator", generator_name)
        
    p_nom = row['capacity_mw']

    if index[0]=='Hydro':
        p_max_pu = (row['estimated_generation_gwh_2017'] * 1000) / (row['capacity_mw'] * 365 * 24 )
    else:
        p_max_pu = 1
        
    n.add("Generator",                           # Adding generators
          generator_name, 
          bus=index[1],
          p_nom_extendable = False,              # The Generators are not extendable
          capital_cost = 0,                      # The Generators have no capital costs
          p_nom=p_nom,
          marginal_cost = marginal_costs[index[0]],
          p_max_pu = p_max_pu,
          carrier=index[0],
         )
    
    for generator_name in n.generators.index:                         # Disregarding existing wind and solar capacities
        if "Wind" in generator_name or "Solar" in generator_name:
            n.remove("Generator", generator_name)
        
display(n.generators.carrier)
Generator
Coal_Košický                Coal
Coal_Trenčiansky            Coal
Gas_Bratislavský             Gas
Gas_Nitriansky               Gas
Gas_Trnavský                 Gas
Hydro_Banskobystrický      Hydro
Hydro_Bratislavský         Hydro
Hydro_Košický              Hydro
Hydro_Trenčiansky          Hydro
Hydro_Trnavský             Hydro
Hydro_Žilinský             Hydro
Nuclear_Nitriansky       Nuclear
Nuclear_Trnavský         Nuclear
Name: carrier, dtype: object
In [22]:
load_time_series = pd.read_csv('load.csv', index_col=0, parse_dates=True)   # Reading the load file
load_data_series_svk = load_time_series[['SK']].copy()                      # Cutting up the Slovakia part of the time series
#display(load_data_series_svk.head(10))

new_index = load_data_series_svk.index.strftime('%Y-%m-%d %H:%M:%S').str.replace('2013', '2017')
load_data_series_svk.index = pd.to_datetime(new_index)
load_data_series_svk = load_data_series_svk.rename_axis('Time')


#display(load_data_series_svk)

# Defining a dictionary including regions population [ref.: https://www.citypopulation.de/en/slovakia/cities/]:
cities_population = {'Banskobystrický':617777,
                     'Bratislavský' : 728370,
                     'Košický': 779505,
                     'Nitriansky': 670696,
                     'Prešovský': 808090,
                     'Trenčiansky': 570675,
                     'Trnavský': 565573,
                     'Žilinský': 688106}

total_population_skv = sum(cities_population.values())    # Calculating the total population of Slovokia

# Adding extra columns as region's load value time series to the primary load time series:
for region in cities_population.keys():
    load_data_series_svk ['skv_load_per_capita'] = load_data_series_svk['SK'] / total_population_skv    # calculating the load per capita
    load_data_series_svk [f'load_{region}'] = load_data_series_svk ['skv_load_per_capita'] * cities_population [region]   # Calculating the load on each region based on the region's population

#display(load_data_series_svk)


for index , bus in n.buses.iterrows():
    Load_name = f'Demand_load_{index}'
    if (n.loads.index == Load_name).any():
        n.remove("Load", Load_name)
    p_set_data = load_data_series_svk[f'load_{index}'] 

    n.add("Load",                           # Adding Electricity loads for each bus (region)
        Load_name, 
        bus=index,
        p_set=p_set_data,             
        carrier="electricity",                     
        )

display(n.loads_t.p_set)

ax = n.loads_t.p_set.iloc[:24].plot(figsize=(10, 6),ylim=[0, 800], ylabel="Load [MW]", xlabel="Time", title='Electricity Demand Load during a single day in different regions of Slovokia (01.01.2017)')
plt.savefig('electricity_demand_load.png')
plt.show()
#n.snapshots
Load Demand_load_Banskobystrický Demand_load_Bratislavský Demand_load_Košický Demand_load_Nitriansky Demand_load_Prešovský Demand_load_Trenčiansky Demand_load_Trnavský Demand_load_Žilinský
snapshot
2017-01-01 00:00:00 364.359299 429.586052 459.745013 395.570447 476.604188 336.578964 333.569851 405.838708
2017-01-01 01:00:00 363.843475 428.977887 459.094152 395.010438 475.929459 336.102469 333.097616 405.264162
2017-01-01 02:00:00 362.529793 427.429033 457.436561 393.584226 474.211084 334.888948 331.894944 403.800928
2017-01-01 03:00:00 365.231299 430.614155 460.845294 396.517143 477.744817 337.384479 334.368164 406.809978
2017-01-01 04:00:00 381.691731 450.021296 481.614908 414.387582 499.276067 352.589897 349.437641 425.144300
... ... ... ... ... ... ... ... ...
2017-12-31 19:00:00 460.693214 543.165440 581.298209 500.156360 602.614826 425.567964 421.763263 513.139474
2017-12-31 20:00:00 435.896276 513.929412 550.009674 473.235307 570.178918 402.661652 399.061740 485.519601
2017-12-31 21:00:00 405.975189 478.651922 512.255539 440.751170 531.040312 375.021879 371.669074 452.192237
2017-12-31 22:00:00 377.655121 445.262061 476.521552 410.005194 493.995935 348.861055 345.742136 420.648154
2017-12-31 23:00:00 372.171151 438.796364 469.601932 404.051465 486.822567 343.795207 340.721578 414.539877

8760 rows × 8 columns

No description has been provided for this image
In [23]:
# Defining regions area as a dictionary[km^2]:
Regions_area = {
                'Banskobystrický' :  9455,
                'Bratislavský' : 2053,
                'Košický': 6753,
                'Nitriansky': 6343,
                'Prešovský' : 8974 ,
                'Trenčiansky': 4502 ,
                'Trnavský': 4145,
                'Žilinský': 6809,
}

# Adding solar and on-shore wind generators per region:
n.buses

for index , bus in n.buses.iterrows():
    solar_generator_name = f'Solar_{index}'
    if solar_generator_name in n.generators.index:
        n.remove("Generator", solar_generator_name)   
    n.add("Generator",                           # Adding solar generators
        solar_generator_name, 
        bus=index,
        p_nom_extendable = False,              # The Generators are not extendable
        capital_cost = 0,                      # The Generators have no capital costs
        marginal_cost = marginalcost_Solar,
        p_nom= Regions_area[index] * cap_per_sqkm,
        p_max_pu = svk_pv_2017_potential_time_series[index],
        carrier='Solar')

    wind_generator_name = f'Wind_{index}'

    if wind_generator_name in n.generators.index:
        n.remove("Generator", wind_generator_name)
    n.add("Generator",                           # Adding wind generators
        wind_generator_name, 
        bus=index,
        p_nom_extendable = False,              # The Generators are not extendable
        capital_cost = 0,                      # The Generators have no capital costs
        marginal_cost = marginalcost_Wind,
        p_nom= Regions_area[index] * cap_per_sqkm,
        p_max_pu = svk_wind_onshore_2017_potential_time_series[index],
        carrier='Wind',
         )


display( n.generators.carrier.map(n.carriers.co2_emissions))
display(n.generators_t.p_max_pu)
Generator
Coal_Košický             0.34
Coal_Trenčiansky         0.34
Gas_Bratislavský         0.20
Gas_Nitriansky           0.20
Gas_Trnavský             0.20
Hydro_Banskobystrický    0.00
Hydro_Bratislavský       0.00
Hydro_Košický            0.00
Hydro_Trenčiansky        0.00
Hydro_Trnavský           0.00
Hydro_Žilinský           0.00
Nuclear_Nitriansky       0.00
Nuclear_Trnavský         0.00
Solar_Banskobystrický    0.00
Wind_Banskobystrický     0.00
Solar_Bratislavský       0.00
Wind_Bratislavský        0.00
Solar_Košický            0.00
Wind_Košický             0.00
Solar_Nitriansky         0.00
Wind_Nitriansky          0.00
Solar_Prešovský          0.00
Wind_Prešovský           0.00
Solar_Trenčiansky        0.00
Wind_Trenčiansky         0.00
Solar_Trnavský           0.00
Wind_Trnavský            0.00
Solar_Žilinský           0.00
Wind_Žilinský            0.00
Name: carrier, dtype: float64
Generator Solar_Banskobystrický Wind_Banskobystrický Solar_Bratislavský Wind_Bratislavský Solar_Košický Wind_Košický Solar_Nitriansky Wind_Nitriansky Solar_Prešovský Wind_Prešovský Solar_Trenčiansky Wind_Trenčiansky Solar_Trnavský Wind_Trnavský Solar_Žilinský Wind_Žilinský
snapshot
2017-01-01 00:00:00 0.0 0.000003 0.0 0.000060 0.0 0.007087 0.0 0.000000e+00 0.0 0.030756 0.0 0.000335 0.0 0.000001 0.0 0.013306
2017-01-01 01:00:00 0.0 0.000000 0.0 0.000005 0.0 0.007612 0.0 0.000000e+00 0.0 0.034854 0.0 0.000475 0.0 0.000040 0.0 0.014443
2017-01-01 02:00:00 0.0 0.000000 0.0 0.000066 0.0 0.007142 0.0 1.130307e-06 0.0 0.039974 0.0 0.000691 0.0 0.000331 0.0 0.015357
2017-01-01 03:00:00 0.0 0.000000 0.0 0.000171 0.0 0.003622 0.0 5.924064e-07 0.0 0.041654 0.0 0.001257 0.0 0.000336 0.0 0.017744
2017-01-01 04:00:00 0.0 0.000000 0.0 0.000144 0.0 0.003500 0.0 0.000000e+00 0.0 0.047368 0.0 0.001554 0.0 0.000341 0.0 0.020749
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
2017-12-31 19:00:00 0.0 0.009669 0.0 0.135380 0.0 0.021920 0.0 9.197736e-02 0.0 0.127549 0.0 0.094018 0.0 0.114936 0.0 0.116168
2017-12-31 20:00:00 0.0 0.006298 0.0 0.135983 0.0 0.022566 0.0 5.302415e-02 0.0 0.130486 0.0 0.072025 0.0 0.096996 0.0 0.096775
2017-12-31 21:00:00 0.0 0.002918 0.0 0.151189 0.0 0.019083 0.0 5.068621e-02 0.0 0.128027 0.0 0.053257 0.0 0.105005 0.0 0.081824
2017-12-31 22:00:00 0.0 0.000884 0.0 0.114648 0.0 0.016145 0.0 4.555578e-02 0.0 0.118321 0.0 0.033475 0.0 0.077719 0.0 0.082753
2017-12-31 23:00:00 0.0 0.000821 0.0 0.158193 0.0 0.023487 0.0 4.725993e-02 0.0 0.119073 0.0 0.027659 0.0 0.121956 0.0 0.073361

8760 rows × 16 columns

In [24]:
# Comparing solar and wind potential in January 2017:

n.generators_t.p_max_pu.groupby(n.generators.carrier, axis=1).mean().iloc[:744].plot(ylabel="p.u.", figsize=(15, 6), color=['y', 'b'])
plt.savefig('Mean renewable energy capacity of implemented Generators-Jan-2017.png')
plt.show()

# Comparing solar and wind potential in 1.January 2017 (as a single day):
n.generators_t.p_max_pu.groupby(n.generators.carrier, axis=1).mean().iloc[:24].plot(ylabel="p.u.", figsize=(15, 6), color=['y', 'b'])
plt.savefig('Mean renewable energy capacity of implemented Generators-one day-1.Jan-2017.png')
plt.show()
No description has been provided for this image
No description has been provided for this image
In [25]:
# Difining a list of tuples, including neighboring regions to create transmission lines between them:
neighbor_buses = [('Bratislavský','Trnavský'),
                  ('Trnavský','Nitriansky'),
                  ('Trnavský','Trenčiansky'),
                  ('Nitriansky','Trenčiansky'),
                  ('Nitriansky','Banskobystrický'),
                  ('Banskobystrický','Trenčiansky'),
                  ('Žilinský','Trenčiansky'),
                  ('Banskobystrický','Žilinský'),
                  ('Prešovský','Žilinský'),
                  ('Prešovský','Banskobystrický'),
                  ('Košický','Banskobystrický'),
                  ('Košický','Prešovský')]

for region1, region2 in neighbor_buses:  
    
    distance = np.sqrt((regions_points[region1].x - regions_points[region2].x) ** 2 +    # Calculating the distance for each transmission line based on the representative points of buses
                       (regions_points[region1].y - regions_points[region2].y) ** 2)   
    length = 1.5 * distance   # As requested, the length of 1.5 times the crow-fly distance between the regions’ representative points
    cost = 400 * length       # Assumption for costs of 400€/MW/km

    # Adding bidirectional links between neighboring regions
    TL1_name = f"TL1_{region1}_{region2}"
    if TL1_name in n.links.index:
        n.remove("Link", TL1_name)
    n.add("Link",
          TL1_name,
          bus0=region1,
          bus1=region2,
          capital_cost=cost,
          p_max_pu=1.0,         # no transmission losses
          efficiency=1.0,       # no transmission losses
          length=length,
         )

    TL2_name = f"TL2_{region2}_{region1}"
    if TL2_name in n.links.index:
        n.remove("Link", TL2_name)
    n.add("Link",
          TL2_name,
          bus0=region2,
          bus1=region1,
          capital_cost=cost,
          p_max_pu=1.0,           # no transmission losses
          efficiency=1.0,         # no transmission losses
          length=length,
         )

display(n.links)
attribute bus0 bus1 type carrier efficiency build_year lifetime p_nom p_nom_extendable p_nom_min ... shut_down_cost min_up_time min_down_time up_time_before down_time_before ramp_limit_up ramp_limit_down ramp_limit_start_up ramp_limit_shut_down p_nom_opt
Link
TL1_Bratislavský_Trnavský Bratislavský Trnavský 1.0 0 inf 0.0 False 0.0 ... 0.0 0 0 1 0 NaN NaN 1.0 1.0 0.0
TL2_Trnavský_Bratislavský Trnavský Bratislavský 1.0 0 inf 0.0 False 0.0 ... 0.0 0 0 1 0 NaN NaN 1.0 1.0 0.0
TL1_Trnavský_Nitriansky Trnavský Nitriansky 1.0 0 inf 0.0 False 0.0 ... 0.0 0 0 1 0 NaN NaN 1.0 1.0 0.0
TL2_Nitriansky_Trnavský Nitriansky Trnavský 1.0 0 inf 0.0 False 0.0 ... 0.0 0 0 1 0 NaN NaN 1.0 1.0 0.0
TL1_Trnavský_Trenčiansky Trnavský Trenčiansky 1.0 0 inf 0.0 False 0.0 ... 0.0 0 0 1 0 NaN NaN 1.0 1.0 0.0
TL2_Trenčiansky_Trnavský Trenčiansky Trnavský 1.0 0 inf 0.0 False 0.0 ... 0.0 0 0 1 0 NaN NaN 1.0 1.0 0.0
TL1_Nitriansky_Trenčiansky Nitriansky Trenčiansky 1.0 0 inf 0.0 False 0.0 ... 0.0 0 0 1 0 NaN NaN 1.0 1.0 0.0
TL2_Trenčiansky_Nitriansky Trenčiansky Nitriansky 1.0 0 inf 0.0 False 0.0 ... 0.0 0 0 1 0 NaN NaN 1.0 1.0 0.0
TL1_Nitriansky_Banskobystrický Nitriansky Banskobystrický 1.0 0 inf 0.0 False 0.0 ... 0.0 0 0 1 0 NaN NaN 1.0 1.0 0.0
TL2_Banskobystrický_Nitriansky Banskobystrický Nitriansky 1.0 0 inf 0.0 False 0.0 ... 0.0 0 0 1 0 NaN NaN 1.0 1.0 0.0
TL1_Banskobystrický_Trenčiansky Banskobystrický Trenčiansky 1.0 0 inf 0.0 False 0.0 ... 0.0 0 0 1 0 NaN NaN 1.0 1.0 0.0
TL2_Trenčiansky_Banskobystrický Trenčiansky Banskobystrický 1.0 0 inf 0.0 False 0.0 ... 0.0 0 0 1 0 NaN NaN 1.0 1.0 0.0
TL1_Žilinský_Trenčiansky Žilinský Trenčiansky 1.0 0 inf 0.0 False 0.0 ... 0.0 0 0 1 0 NaN NaN 1.0 1.0 0.0
TL2_Trenčiansky_Žilinský Trenčiansky Žilinský 1.0 0 inf 0.0 False 0.0 ... 0.0 0 0 1 0 NaN NaN 1.0 1.0 0.0
TL1_Banskobystrický_Žilinský Banskobystrický Žilinský 1.0 0 inf 0.0 False 0.0 ... 0.0 0 0 1 0 NaN NaN 1.0 1.0 0.0
TL2_Žilinský_Banskobystrický Žilinský Banskobystrický 1.0 0 inf 0.0 False 0.0 ... 0.0 0 0 1 0 NaN NaN 1.0 1.0 0.0
TL1_Prešovský_Žilinský Prešovský Žilinský 1.0 0 inf 0.0 False 0.0 ... 0.0 0 0 1 0 NaN NaN 1.0 1.0 0.0
TL2_Žilinský_Prešovský Žilinský Prešovský 1.0 0 inf 0.0 False 0.0 ... 0.0 0 0 1 0 NaN NaN 1.0 1.0 0.0
TL1_Prešovský_Banskobystrický Prešovský Banskobystrický 1.0 0 inf 0.0 False 0.0 ... 0.0 0 0 1 0 NaN NaN 1.0 1.0 0.0
TL2_Banskobystrický_Prešovský Banskobystrický Prešovský 1.0 0 inf 0.0 False 0.0 ... 0.0 0 0 1 0 NaN NaN 1.0 1.0 0.0
TL1_Košický_Banskobystrický Košický Banskobystrický 1.0 0 inf 0.0 False 0.0 ... 0.0 0 0 1 0 NaN NaN 1.0 1.0 0.0
TL2_Banskobystrický_Košický Banskobystrický Košický 1.0 0 inf 0.0 False 0.0 ... 0.0 0 0 1 0 NaN NaN 1.0 1.0 0.0
TL1_Košický_Prešovský Košický Prešovský 1.0 0 inf 0.0 False 0.0 ... 0.0 0 0 1 0 NaN NaN 1.0 1.0 0.0
TL2_Prešovský_Košický Prešovský Košický 1.0 0 inf 0.0 False 0.0 ... 0.0 0 0 1 0 NaN NaN 1.0 1.0 0.0

24 rows × 32 columns

In [26]:
# Adding Battery storage units:

for bus in n.buses.index:
    battery_storage_U_name = f"Battery_storage_{bus}"
    if battery_storage_U_name in n.storage_units.index:
        n.remove("StorageUnit", battery_storage_U_name)
    n.add(
        "StorageUnit",
        battery_storage_U_name,
        max_hours= 6,
        capital_cost = 6 * capitalcost_battery_storage + capitalcost_battery_inverter,
        bus = bus,
        marginal_cost =0,
        carrier = "battery storage",
        efficiency_dispatch = df_battery_inverter.loc['efficiency', 'value'],
        efficiency_store = df_battery_inverter.loc['efficiency', 'value'],
        cyclic_state_of_charge = True,
        p_nom_extendable = True,
  
    )

n.storage_units
Out[26]:
attribute bus control type p_nom p_nom_extendable p_nom_min p_nom_max p_min_pu p_max_pu p_set ... state_of_charge_initial_per_period state_of_charge_set cyclic_state_of_charge cyclic_state_of_charge_per_period max_hours efficiency_store efficiency_dispatch standing_loss inflow p_nom_opt
StorageUnit
Battery_storage_Banskobystrický Banskobystrický PQ 0.0 True 0.0 inf -1.0 1.0 0.0 ... False NaN True True 6.0 0.96 0.96 0.0 0.0 0.0
Battery_storage_Bratislavský Bratislavský PQ 0.0 True 0.0 inf -1.0 1.0 0.0 ... False NaN True True 6.0 0.96 0.96 0.0 0.0 0.0
Battery_storage_Košický Košický PQ 0.0 True 0.0 inf -1.0 1.0 0.0 ... False NaN True True 6.0 0.96 0.96 0.0 0.0 0.0
Battery_storage_Nitriansky Nitriansky PQ 0.0 True 0.0 inf -1.0 1.0 0.0 ... False NaN True True 6.0 0.96 0.96 0.0 0.0 0.0
Battery_storage_Prešovský Prešovský PQ 0.0 True 0.0 inf -1.0 1.0 0.0 ... False NaN True True 6.0 0.96 0.96 0.0 0.0 0.0
Battery_storage_Trenčiansky Trenčiansky PQ 0.0 True 0.0 inf -1.0 1.0 0.0 ... False NaN True True 6.0 0.96 0.96 0.0 0.0 0.0
Battery_storage_Trnavský Trnavský PQ 0.0 True 0.0 inf -1.0 1.0 0.0 ... False NaN True True 6.0 0.96 0.96 0.0 0.0 0.0
Battery_storage_Žilinský Žilinský PQ 0.0 True 0.0 inf -1.0 1.0 0.0 ... False NaN True True 6.0 0.96 0.96 0.0 0.0 0.0

8 rows × 29 columns

In [27]:
# Adding hydrogen storage units

for bus in n.buses.index:
    storage_U_name = f"hydrogen_storage_{bus}"
    if storage_U_name in n.storage_units.index:
        n.remove("StorageUnit", storage_U_name)
    n.add(
        "StorageUnit",
        storage_U_name,
        max_hours= 336,
        capital_cost = 336 * capitalcost_hydrogen_storage_underground+capitalcost_fuel_cell+capitalcost_electrolysis,
        bus = bus,
        marginal_cost = 0,
        carrier = "hydrogen storage underground",
        efficiency_dispatch = df_fuel_cell.loc['efficiency', 'value'],
        efficiency_store = df_electrolysis.loc['efficiency', 'value'],
        cyclic_state_of_charge = True,
        p_nom_extendable = True,   
    )

display(n.storage_units)
attribute bus control type p_nom p_nom_extendable p_nom_min p_nom_max p_min_pu p_max_pu p_set ... state_of_charge_initial_per_period state_of_charge_set cyclic_state_of_charge cyclic_state_of_charge_per_period max_hours efficiency_store efficiency_dispatch standing_loss inflow p_nom_opt
StorageUnit
Battery_storage_Banskobystrický Banskobystrický PQ 0.0 True 0.0 inf -1.0 1.0 0.0 ... False NaN True True 6.0 0.9600 0.96 0.0 0.0 0.0
Battery_storage_Bratislavský Bratislavský PQ 0.0 True 0.0 inf -1.0 1.0 0.0 ... False NaN True True 6.0 0.9600 0.96 0.0 0.0 0.0
Battery_storage_Košický Košický PQ 0.0 True 0.0 inf -1.0 1.0 0.0 ... False NaN True True 6.0 0.9600 0.96 0.0 0.0 0.0
Battery_storage_Nitriansky Nitriansky PQ 0.0 True 0.0 inf -1.0 1.0 0.0 ... False NaN True True 6.0 0.9600 0.96 0.0 0.0 0.0
Battery_storage_Prešovský Prešovský PQ 0.0 True 0.0 inf -1.0 1.0 0.0 ... False NaN True True 6.0 0.9600 0.96 0.0 0.0 0.0
Battery_storage_Trenčiansky Trenčiansky PQ 0.0 True 0.0 inf -1.0 1.0 0.0 ... False NaN True True 6.0 0.9600 0.96 0.0 0.0 0.0
Battery_storage_Trnavský Trnavský PQ 0.0 True 0.0 inf -1.0 1.0 0.0 ... False NaN True True 6.0 0.9600 0.96 0.0 0.0 0.0
Battery_storage_Žilinský Žilinský PQ 0.0 True 0.0 inf -1.0 1.0 0.0 ... False NaN True True 6.0 0.9600 0.96 0.0 0.0 0.0
hydrogen_storage_Banskobystrický Banskobystrický PQ 0.0 True 0.0 inf -1.0 1.0 0.0 ... False NaN True True 336.0 0.6532 0.50 0.0 0.0 0.0
hydrogen_storage_Bratislavský Bratislavský PQ 0.0 True 0.0 inf -1.0 1.0 0.0 ... False NaN True True 336.0 0.6532 0.50 0.0 0.0 0.0
hydrogen_storage_Košický Košický PQ 0.0 True 0.0 inf -1.0 1.0 0.0 ... False NaN True True 336.0 0.6532 0.50 0.0 0.0 0.0
hydrogen_storage_Nitriansky Nitriansky PQ 0.0 True 0.0 inf -1.0 1.0 0.0 ... False NaN True True 336.0 0.6532 0.50 0.0 0.0 0.0
hydrogen_storage_Prešovský Prešovský PQ 0.0 True 0.0 inf -1.0 1.0 0.0 ... False NaN True True 336.0 0.6532 0.50 0.0 0.0 0.0
hydrogen_storage_Trenčiansky Trenčiansky PQ 0.0 True 0.0 inf -1.0 1.0 0.0 ... False NaN True True 336.0 0.6532 0.50 0.0 0.0 0.0
hydrogen_storage_Trnavský Trnavský PQ 0.0 True 0.0 inf -1.0 1.0 0.0 ... False NaN True True 336.0 0.6532 0.50 0.0 0.0 0.0
hydrogen_storage_Žilinský Žilinský PQ 0.0 True 0.0 inf -1.0 1.0 0.0 ... False NaN True True 336.0 0.6532 0.50 0.0 0.0 0.0

16 rows × 29 columns

In [28]:
# Adding shedding generator:

n.madd(
            "Generator",
            n.buses.index,
            suffix=" load shedding",
            bus=n.buses.index,
            carrier="load",
            sign=1e-3,  # Adjust sign to measure p and p_nom in kW instead of MW
            marginal_cost=1000,  # Eur/kWh
            p_nom=1e9,  # kW
)
Out[28]:
Index(['Banskobystrický load shedding', 'Bratislavský load shedding',
       'Košický load shedding', 'Nitriansky load shedding',
       'Prešovský load shedding', 'Trenčiansky load shedding',
       'Trnavský load shedding', 'Žilinský load shedding'],
      dtype='object', name='Bus')

1) Solving The Model without a limit on CO2 emissions:¶

In [29]:
# Calculating the CO2 emissions for each region based on the resulted power output:

e = (
    n.generators_t.p
    / n.generators.efficiency
    * n.generators.carrier.map(n.carriers.co2_emissions)
)
e.tail(20)

# Calculating the CO2 emissions in total (for the whole country)[tonnes/MWh thermal]:
display(e.sum().sum())

# Adding global constraint for CO2 emission 
n.add(
    "GlobalConstraint",
    "emission_limit",
    carrier_attribute="co2_emissions",
    sense="<=",
    constant=e.sum().sum() * 1,               # no reduction constraints in CO2 emission 
)


# Solving the Model:

n.optimize(solver_name="highs")

n.generators_t.p    # Retrieveing the power output
0.0
INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.io:Writing objective.
Writing constraints.: 100%|████████████████████████████████████████████████████████████| 15/15 [00:13<00:00,  1.08it/s]
Writing continuous variables.: 100%|█████████████████████████████████████████████████████| 7/7 [00:02<00:00,  2.52it/s]
INFO:linopy.io: Writing time: 17.16s
INFO:linopy.solvers:Log file at C:\Users\Mahtab\AppData\Local\Temp\highs.log.
INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 954856 primals, 2119937 duals
Objective: 4.59e+08
Solver model: available
Solver message: optimal

INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-fix-p-lower, Generator-fix-p-upper, Link-fix-p-lower, Link-fix-p-upper, StorageUnit-ext-p_dispatch-lower, StorageUnit-ext-p_dispatch-upper, StorageUnit-ext-p_store-lower, StorageUnit-ext-p_store-upper, StorageUnit-ext-state_of_charge-lower, StorageUnit-ext-state_of_charge-upper, StorageUnit-energy_balance were not assigned to the network.
Out[29]:
Generator Coal_Košický Coal_Trenčiansky Gas_Bratislavský Gas_Nitriansky Gas_Trnavský Hydro_Banskobystrický Hydro_Bratislavský Hydro_Košický Hydro_Trenčiansky Hydro_Trnavský ... Solar_Žilinský Wind_Žilinský Banskobystrický load shedding Bratislavský load shedding Košický load shedding Nitriansky load shedding Prešovský load shedding Trenčiansky load shedding Trnavský load shedding Žilinský load shedding
snapshot
2017-01-01 00:00:00 -0.0 -0.0 -0.0 -0.0 -0.0 14.381279 14.923516 18.437215 6.369863 247.155251 ... -0.0 271.806237 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0
2017-01-01 01:00:00 -0.0 -0.0 -0.0 -0.0 -0.0 14.381279 14.923516 18.437215 6.369863 247.155251 ... -0.0 295.034377 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0
2017-01-01 02:00:00 -0.0 -0.0 -0.0 -0.0 -0.0 14.381279 14.923516 18.437215 6.369863 247.155251 ... -0.0 313.706895 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0
2017-01-01 03:00:00 -0.0 -0.0 -0.0 -0.0 -0.0 14.381279 14.923516 18.437215 6.369863 247.155251 ... -0.0 362.456009 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0
2017-01-01 04:00:00 -0.0 -0.0 -0.0 -0.0 -0.0 14.381279 14.923516 18.437215 6.369863 247.155251 ... -0.0 423.834183 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
2017-12-31 19:00:00 -0.0 -0.0 -0.0 -0.0 -0.0 14.381279 -0.000000 18.437215 -0.000000 -0.000000 ... -0.0 -0.000000 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0
2017-12-31 20:00:00 -0.0 -0.0 -0.0 -0.0 -0.0 14.381279 14.923516 18.437215 6.369863 -0.000000 ... -0.0 485.519601 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0
2017-12-31 21:00:00 -0.0 -0.0 -0.0 -0.0 -0.0 14.381279 14.923516 18.437215 6.369863 -0.000000 ... -0.0 -0.000000 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0
2017-12-31 22:00:00 -0.0 -0.0 -0.0 -0.0 -0.0 14.381279 14.923516 18.437215 6.369863 -0.000000 ... -0.0 -0.000000 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0
2017-12-31 23:00:00 -0.0 -0.0 -0.0 -0.0 -0.0 14.381279 -0.000000 18.437215 6.369863 -0.000000 ... -0.0 414.539877 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0

8760 rows × 37 columns

In [30]:
# Regional distribution of loads in Slovakia


load = n.loads_t.p_set.sum(axis=0).groupby(n.loads.bus).sum()
display(load)

fig = plt.figure(figsize=(15, 6))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.EqualEarth())

# Plotting the network with bus sizes proportional to the load
n.plot(
    ax=ax,
    bus_sizes=load / 2e9,  
    geomap=True, 
    color_geomap=True,  
    projection=ccrs.EqualEarth()
)
svk_regions.plot(ax=ax, facecolor='lightblue', edgecolor='black', linewidth=1.5, transform=ccrs.PlateCarree())
ax.set_extent([minx, maxx, miny, maxy], crs=ccrs.PlateCarree())  
plt.title('Regional Distribution of Load')
plt.show()
bus
Banskobystrický    3.480420e+06
Bratislavský       4.103476e+06
Košický            4.391560e+06
Nitriansky         3.778554e+06
Prešovský          4.552601e+06
Trenčiansky        3.215057e+06
Trnavský           3.186314e+06
Žilinský           3.876638e+06
dtype: float64
No description has been provided for this image
In [31]:
carrier_to_exclude = 'load'                # Excluding the "load" carriers, which are created during the Shedding Generator defining
load_excluded_generators = n.generators[n.generators.carrier != carrier_to_exclude]

dispatch_by_carrier = n.generators_t.p.groupby(load_excluded_generators.carrier, axis=1).sum().div(1e3)
#display(p_by_carrier.head(20))

fig, ax = plt.subplots(figsize=(15, 8))

dispatch_by_carrier.iloc[:720].plot(
    kind="area",
    ax=ax,
    linewidth=0,
    cmap="tab20b",
)

ax.legend(ncol=5, loc="upper left", frameon=False)

ax.set_ylabel("GW")
ax.set_title('Hourly dispatch per technology- JAN 2017')
ax.set_ylim(0, 9);

#########################################################
fig, ax = plt.subplots(figsize=(15, 8))

dispatch_by_carrier.iloc[:24].plot(
    kind="area",
    ax=ax,
    linewidth=0,
    cmap="tab20b",
)

ax.legend(ncol=5, loc="upper left", frameon=False)

ax.set_ylabel("GW")
ax.set_title('Hourly dispatch per technology through a single day- 1.JAN 2017')

ax.set_ylim(0, 9);
No description has been provided for this image
No description has been provided for this image
In [32]:
# The aggregate dispatch of the pumped hydro storage units and the state of charge throughout the day:

fig, ax = plt.subplots(figsize=(15, 8))

battery_storage_units = n.storage_units[n.storage_units.carrier == 'battery storage']

p_battery_storage = n.storage_units_t.p[battery_storage_units.index].sum(axis=1).div(1e3)
state_of_charge_battery = n.storage_units_t.state_of_charge[battery_storage_units.index].sum(axis=1).div(1e3)

hydro_storage_units = n.storage_units[n.storage_units.carrier == 'hydrogen storage underground']
p_hydro_storage = n.storage_units_t.p[hydro_storage_units.index].sum(axis=1).div(1e3)
state_of_charge_hydro = n.storage_units_t.state_of_charge[hydro_storage_units.index].sum(axis=1).div(1e3)

total_storage_region1 = p_hydro_storage + p_battery_storage
total_SOC_region1 = state_of_charge_hydro +  state_of_charge_battery
#display(p_hydro_storage.iloc[:24])
total_storage_region1.iloc[:24].plot(label=" dispatch [GW]", ax=ax)
total_SOC_region1.iloc[:24].plot(label=" state of charge [GWh]", ax=ax)

ax.grid()
ax.legend()
ax.set_ylabel("MWh or MW")
ax.set_title('Total Storage: Dispatch and State of Charge throughout the day in Banskobystrický ')

plt.show()
###############################################################

# The aggregate dispatch of the pumped hydro storage units and the state of charge throughout a Month (January):

fig, ax = plt.subplots(figsize=(15, 8))


total_storage_region1.iloc[:720].plot(label=" dispatch [GW]", ax=ax)
total_SOC_region1.iloc[:720].plot(label=" state of charge [GWh]", ax=ax)

ax.grid()
ax.legend()
ax.set_ylabel("MWh or MW")
ax.set_title('Total Storage: Dispatch and State of Charge throughout a month in Banskobystrický ')

plt.show()



##############################################################
# the aggregate dispatch of the pumped hydro storage units and the state of charge throughout the whole year:

fig, ax = plt.subplots(figsize=(15, 8))


total_storage_region1.plot(label=" dispatch [GW]", ax=ax)
total_SOC_region1.plot(label=" state of charge [GWh]", ax=ax)

ax.grid()
ax.legend()
ax.set_ylabel("MWh or MW")
ax.set_title('Total Storage: Dispatch and State of Charge throughout the whole year in Banskobystrický ')

plt.show()
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
In [33]:
# Locational marginal prices:


fig = plt.figure(figsize=(7, 7))
ax = plt.axes(projection=ccrs.EqualEarth())

norm = plt.Normalize(vmin=0, vmax=50)  # €/MWh

n.plot(
    ax=ax,
    bus_colors=n.buses_t.marginal_price.mean(),
    bus_cmap="plasma",
    bus_norm=norm,

)

cax = fig.add_axes([0.95, 0.15, 0.02, 0.7])  
cb = plt.colorbar(
    plt.cm.ScalarMappable(cmap="plasma", norm=norm),
    label="LMP [€/MWh]",
    shrink=0.6,
    cax=cax,  
)
ax.set_extent([minx, maxx, miny, maxy], crs=ccrs.PlateCarree())  

plt.show()
No description has been provided for this image
In [34]:
display(n.generators['marginal_cost'])
generation_cost = (n.generators_t.p * n.generators['marginal_cost']).sum(axis=1)

sorted_costs = generation_cost[np.argsort(-generation_cost)]    # Sorting generation costs in ascending order
cumulative_probability = np.arange(1, len(sorted_costs) + 1) / len(sorted_costs)   # Calculating the cumulative probability

# Plotting:
plt.figure(figsize=(10, 6))
plt.plot(cumulative_probability * 100, sorted_costs, linestyle='-')
plt.xlabel('Cumulative Probability (%)')
plt.ylabel(' price (€)')
plt.title('Cost Duration Curve')
plt.grid(True)
plt.show()
Generator
Coal_Košický                       31.431739
Coal_Trenčiansky                   31.431739
Gas_Bratislavský                   44.520508
Gas_Nitriansky                     44.520508
Gas_Trnavský                       44.520508
Hydro_Banskobystrický               0.000000
Hydro_Bratislavský                  0.000000
Hydro_Košický                       0.000000
Hydro_Trenčiansky                   0.000000
Hydro_Trnavský                      0.000000
Hydro_Žilinský                      0.000000
Nuclear_Nitriansky                 13.778923
Nuclear_Trnavský                   13.778923
Solar_Banskobystrický               0.010000
Wind_Banskobystrický                0.000000
Solar_Bratislavský                  0.010000
Wind_Bratislavský                   0.000000
Solar_Košický                       0.010000
Wind_Košický                        0.000000
Solar_Nitriansky                    0.010000
Wind_Nitriansky                     0.000000
Solar_Prešovský                     0.010000
Wind_Prešovský                      0.000000
Solar_Trenčiansky                   0.010000
Wind_Trenčiansky                    0.000000
Solar_Trnavský                      0.010000
Wind_Trnavský                       0.000000
Solar_Žilinský                      0.010000
Wind_Žilinský                       0.000000
Banskobystrický load shedding    1000.000000
Bratislavský load shedding       1000.000000
Košický load shedding            1000.000000
Nitriansky load shedding         1000.000000
Prešovský load shedding          1000.000000
Trenčiansky load shedding        1000.000000
Trnavský load shedding           1000.000000
Žilinský load shedding           1000.000000
Name: marginal_cost, dtype: float64
No description has been provided for this image
In [35]:
n.global_constraints.mu
Out[35]:
GlobalConstraint
emission_limit   -28014.414082
Name: mu, dtype: float64

2) Solving the model with a CO2 emission reduction of 100% (no emissions):¶

In [36]:
#Setting a constraint to  solver:
n.global_constraints.loc["emission_limit", "constant"] = 0
In [37]:
 n.optimize(solver_name="highs")
INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.io:Writing objective.
Writing constraints.: 100%|████████████████████████████████████████████████████████████| 15/15 [00:15<00:00,  1.06s/it]
Writing continuous variables.: 100%|█████████████████████████████████████████████████████| 7/7 [00:02<00:00,  2.76it/s]
INFO:linopy.io: Writing time: 18.9s
INFO:linopy.solvers:Log file at C:\Users\Mahtab\AppData\Local\Temp\highs.log.
INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 954856 primals, 2119937 duals
Objective: 4.59e+08
Solver model: available
Solver message: optimal

INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-fix-p-lower, Generator-fix-p-upper, Link-fix-p-lower, Link-fix-p-upper, StorageUnit-ext-p_dispatch-lower, StorageUnit-ext-p_dispatch-upper, StorageUnit-ext-p_store-lower, StorageUnit-ext-p_store-upper, StorageUnit-ext-state_of_charge-lower, StorageUnit-ext-state_of_charge-upper, StorageUnit-energy_balance were not assigned to the network.
Out[37]:
('ok', 'optimal')
In [38]:
n.generators_t.p
Out[38]:
Generator Coal_Košický Coal_Trenčiansky Gas_Bratislavský Gas_Nitriansky Gas_Trnavský Hydro_Banskobystrický Hydro_Bratislavský Hydro_Košický Hydro_Trenčiansky Hydro_Trnavský ... Solar_Žilinský Wind_Žilinský Banskobystrický load shedding Bratislavský load shedding Košický load shedding Nitriansky load shedding Prešovský load shedding Trenčiansky load shedding Trnavský load shedding Žilinský load shedding
snapshot
2017-01-01 00:00:00 -0.0 -0.0 -0.0 -0.0 -0.0 14.381279 14.923516 18.437215 6.369863 247.155251 ... -0.0 271.806237 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0
2017-01-01 01:00:00 -0.0 -0.0 -0.0 -0.0 -0.0 14.381279 14.923516 18.437215 6.369863 247.155251 ... -0.0 295.034377 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0
2017-01-01 02:00:00 -0.0 -0.0 -0.0 -0.0 -0.0 14.381279 14.923516 18.437215 6.369863 247.155251 ... -0.0 313.706895 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0
2017-01-01 03:00:00 -0.0 -0.0 -0.0 -0.0 -0.0 14.381279 14.923516 18.437215 6.369863 247.155251 ... -0.0 362.456009 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0
2017-01-01 04:00:00 -0.0 -0.0 -0.0 -0.0 -0.0 14.381279 14.923516 18.437215 6.369863 247.155251 ... -0.0 423.834183 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
2017-12-31 19:00:00 -0.0 -0.0 -0.0 -0.0 -0.0 14.381279 -0.000000 18.437215 -0.000000 -0.000000 ... -0.0 -0.000000 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0
2017-12-31 20:00:00 -0.0 -0.0 -0.0 -0.0 -0.0 14.381279 14.923516 18.437215 6.369863 -0.000000 ... -0.0 485.519601 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0
2017-12-31 21:00:00 -0.0 -0.0 -0.0 -0.0 -0.0 14.381279 14.923516 18.437215 6.369863 -0.000000 ... -0.0 -0.000000 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0
2017-12-31 22:00:00 -0.0 -0.0 -0.0 -0.0 -0.0 14.381279 14.923516 18.437215 6.369863 -0.000000 ... -0.0 -0.000000 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0
2017-12-31 23:00:00 -0.0 -0.0 -0.0 -0.0 -0.0 14.381279 -0.000000 18.437215 6.369863 -0.000000 ... -0.0 414.539877 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0

8760 rows × 37 columns

In [39]:
carrier_to_exclude = 'load'                # Excluding the "load" carriers, which are created during the Shedding Generator defining
load_excluded_generators = n.generators[n.generators.carrier != carrier_to_exclude]

dispatch_by_carrier = n.generators_t.p.groupby(load_excluded_generators.carrier, axis=1).sum().div(1e3)
#display(p_by_carrier.head(20))

fig, ax = plt.subplots(figsize=(15, 8))

dispatch_by_carrier.iloc[:720].plot(
    kind="area",
    ax=ax,
    linewidth=0,
    cmap="tab20b",
)

ax.legend(ncol=5, loc="upper left", frameon=False)

ax.set_ylabel("GW")
ax.set_title('Hourly dispatch per technology- JAN 2017')
ax.set_ylim(0, 9);

#########################################################
fig, ax = plt.subplots(figsize=(15, 8))

dispatch_by_carrier.iloc[:24].plot(
    kind="area",
    ax=ax,
    linewidth=0,
    cmap="tab20b",
)

ax.legend(ncol=5, loc="upper left", frameon=False)

ax.set_ylabel("GW")
ax.set_title('Hourly dispatch per technology through a single day- 1.JAN 2017')

ax.set_ylim(0, 9);
No description has been provided for this image
No description has been provided for this image
In [40]:
generation_cost = (n.generators_t.p * n.generators['marginal_cost']).sum(axis=1)

sorted_costs = generation_cost[np.argsort(-generation_cost)]   # Sorting generation costs in ascending order

cumulative_probability = np.arange(1, len(sorted_costs) + 1) / len(sorted_costs)   # Calculating the cumulative probability


# Ploting:
plt.figure(figsize=(10, 6))
plt.plot(cumulative_probability * 100, sorted_costs, linestyle='-')
plt.xlabel('Cumulative Probability (%)')
plt.ylabel(' price (€)')
plt.title('Cost Duration Curve')
plt.grid(True)
plt.show()
No description has been provided for this image
In [41]:
n.global_constraints.mu
Out[41]:
GlobalConstraint
emission_limit   -28014.414082
Name: mu, dtype: float64
In [42]:
# Locational marginal prices:

fig = plt.figure(figsize=(7, 7))
ax = plt.axes(projection=ccrs.EqualEarth())

norm = plt.Normalize(vmin=0, vmax=50)  # €/MWh

n.plot(
    ax=ax,
    bus_colors=n.buses_t.marginal_price.mean(),
    bus_cmap="plasma",
    bus_norm=norm,

)

cax = fig.add_axes([0.95, 0.15, 0.02, 0.7])  
cb = plt.colorbar(
    plt.cm.ScalarMappable(cmap="plasma", norm=norm),
    label="LMP [€/MWh]",
    shrink=0.6,
    cax=cax,  
)
ax.set_extent([minx, maxx, miny, maxy], crs=ccrs.PlateCarree())  

plt.show()
No description has been provided for this image
In [43]:
# Reduction wind and solar potential in 25% steps:

for reduction_step in np.arange(0.75, -0.01, -0.25):
    
    previou_step = reduction_step + 0.25
    reduction_factor =  reduction_step / previou_step
    
    for generator_name, generator_data in n.generators.iterrows():
        if generator_data['carrier'] in ['Solar', 'Wind']:
           n.generators.at[generator_name, 'p_nom'] *= reduction_factor
        
    n.optimize(solver_name="highs")
    print(f'{(1- round(reduction_step * 100))} % Reduction of Renewable Energy installable potential')
    n.optimize(solver_name="highs")
    print()
    print(f'CO2 Shadow price in case of {(1- round(reduction_step * 100))} % Reduction of Renewable Energy installable potential:')
    print()
    display(n.global_constraints.mu)
    print()
    ###################################################################################################
    # Dispatch power per technology :
    carrier_to_exclude = 'load'                # Excluding the "load" carriers, which are created during the Shedding Generator defining
    load_excluded_generators = n.generators[n.generators.carrier != carrier_to_exclude]
    
    dispatch_by_carrier = n.generators_t.p.groupby(load_excluded_generators.carrier, axis=1).sum().div(1e3)
    
    fig, ax = plt.subplots(figsize=(15, 8))
    
    dispatch_by_carrier.iloc[:720].plot(
        kind="area",
        ax=ax,
        linewidth=0,
        cmap="tab20b",
    )
    
    ax.legend(ncol=5, loc="upper left", frameon=False)
    
    ax.set_ylabel("GW")
    ax.set_title(f'Hourly dispatch per technology ({(1- round(reduction_step * 100))} % Reduction of RE potential)')
    ax.set_ylim(0, 9);
    
    #-----------------------------------------------------------
    fig, ax = plt.subplots(figsize=(15, 8))
    
    dispatch_by_carrier.iloc[:24].plot(
        kind="area",
        ax=ax,
        linewidth=0,
        cmap="tab20b",
    )
    
    ax.legend(ncol=5, loc="upper left", frameon=False)
    
    ax.set_ylabel("GW")
    ax.set_title(f'Hourly dispatch per technology through a single day ({(1- round(reduction_step * 100))} % Reduction of RE potential)')
    
    ax.set_ylim(0, 9);

    #######################################################################################
       
    generation_cost = (n.generators_t.p * n.generators['marginal_cost']).sum(axis=1)
    sorted_costs = generation_cost[np.argsort(-generation_cost)]

    cumulative_probability = np.arange(1, len(sorted_costs) + 1) / len(sorted_costs)
    

    plt.figure(figsize=(10, 6))
    plt.plot(cumulative_probability * 100, sorted_costs, linestyle='-')
    plt.xlabel('Cumulative Probability (%)')
    plt.ylabel(' price (€)')
    plt.title(f'Cost Duration Curve ({(1- round(reduction_step * 100))} % Reduction of RE potential)')
    plt.grid(True)
    plt.show()
    ##############################################################################
    fig = plt.figure(figsize=(7, 7))
    ax = plt.axes(projection=ccrs.EqualEarth())
    
    norm = plt.Normalize(vmin=0, vmax=50)  # €/MWh
    
    n.plot(
        ax=ax,
        bus_colors=n.buses_t.marginal_price.mean(),
        bus_cmap="plasma",
        bus_norm=norm,
    
    )
    
    cax = fig.add_axes([0.95, 0.15, 0.02, 0.7])  
    cb = plt.colorbar(
        plt.cm.ScalarMappable(cmap="plasma", norm=norm),
        label="LMP [€/MWh]",
        shrink=0.6,
        cax=cax,  
    )
    ax.set_extent([minx, maxx, miny, maxy], crs=ccrs.PlateCarree())  
    ax.set_title(f'Electricity price per region ({(1- round(reduction_step * 100))} % Reduction of RE potential)')

    plt.show()
INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.io:Writing objective.
Writing constraints.: 100%|████████████████████████████████████████████████████████████| 15/15 [00:13<00:00,  1.10it/s]
Writing continuous variables.: 100%|█████████████████████████████████████████████████████| 7/7 [00:02<00:00,  2.43it/s]
INFO:linopy.io: Writing time: 17.02s
INFO:linopy.solvers:Log file at C:\Users\Mahtab\AppData\Local\Temp\highs.log.
INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 954856 primals, 2119937 duals
Objective: 4.84e+08
Solver model: available
Solver message: optimal

INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-fix-p-lower, Generator-fix-p-upper, Link-fix-p-lower, Link-fix-p-upper, StorageUnit-ext-p_dispatch-lower, StorageUnit-ext-p_dispatch-upper, StorageUnit-ext-p_store-lower, StorageUnit-ext-p_store-upper, StorageUnit-ext-state_of_charge-lower, StorageUnit-ext-state_of_charge-upper, StorageUnit-energy_balance were not assigned to the network.
-74 % Reduction of Renewable Energy installable potential
INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.io:Writing objective.
Writing constraints.: 100%|████████████████████████████████████████████████████████████| 15/15 [00:19<00:00,  1.27s/it]
Writing continuous variables.: 100%|█████████████████████████████████████████████████████| 7/7 [00:03<00:00,  1.77it/s]
INFO:linopy.io: Writing time: 23.68s
INFO:linopy.solvers:Log file at C:\Users\Mahtab\AppData\Local\Temp\highs.log.
INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 954856 primals, 2119937 duals
Objective: 4.84e+08
Solver model: available
Solver message: optimal

INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-fix-p-lower, Generator-fix-p-upper, Link-fix-p-lower, Link-fix-p-upper, StorageUnit-ext-p_dispatch-lower, StorageUnit-ext-p_dispatch-upper, StorageUnit-ext-p_store-lower, StorageUnit-ext-p_store-upper, StorageUnit-ext-state_of_charge-lower, StorageUnit-ext-state_of_charge-upper, StorageUnit-energy_balance were not assigned to the network.
CO2 Shadow price in case of -74 % Reduction of Renewable Energy installable potential:

GlobalConstraint
emission_limit   -28016.106241
Name: mu, dtype: float64

No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.io:Writing objective.
Writing constraints.: 100%|████████████████████████████████████████████████████████████| 15/15 [00:41<00:00,  2.74s/it]
Writing continuous variables.: 100%|█████████████████████████████████████████████████████| 7/7 [00:08<00:00,  1.17s/it]
INFO:linopy.io: Writing time: 50.74s
INFO:linopy.solvers:Log file at C:\Users\Mahtab\AppData\Local\Temp\highs.log.
INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 954856 primals, 2119937 duals
Objective: 5.32e+08
Solver model: available
Solver message: optimal

INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-fix-p-lower, Generator-fix-p-upper, Link-fix-p-lower, Link-fix-p-upper, StorageUnit-ext-p_dispatch-lower, StorageUnit-ext-p_dispatch-upper, StorageUnit-ext-p_store-lower, StorageUnit-ext-p_store-upper, StorageUnit-ext-state_of_charge-lower, StorageUnit-ext-state_of_charge-upper, StorageUnit-energy_balance were not assigned to the network.
-49 % Reduction of Renewable Energy installable potential
INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.io:Writing objective.
Writing constraints.: 100%|████████████████████████████████████████████████████████████| 15/15 [00:39<00:00,  2.64s/it]
Writing continuous variables.: 100%|█████████████████████████████████████████████████████| 7/7 [00:07<00:00,  1.08s/it]
INFO:linopy.io: Writing time: 48.54s
INFO:linopy.solvers:Log file at C:\Users\Mahtab\AppData\Local\Temp\highs.log.
INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 954856 primals, 2119937 duals
Objective: 5.32e+08
Solver model: available
Solver message: optimal

INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-fix-p-lower, Generator-fix-p-upper, Link-fix-p-lower, Link-fix-p-upper, StorageUnit-ext-p_dispatch-lower, StorageUnit-ext-p_dispatch-upper, StorageUnit-ext-p_store-lower, StorageUnit-ext-p_store-upper, StorageUnit-ext-state_of_charge-lower, StorageUnit-ext-state_of_charge-upper, StorageUnit-energy_balance were not assigned to the network.
CO2 Shadow price in case of -49 % Reduction of Renewable Energy installable potential:

GlobalConstraint
emission_limit   -19881.836628
Name: mu, dtype: float64

No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.io:Writing objective.
Writing constraints.: 100%|████████████████████████████████████████████████████████████| 15/15 [00:39<00:00,  2.65s/it]
Writing continuous variables.: 100%|█████████████████████████████████████████████████████| 7/7 [00:08<00:00,  1.15s/it]
INFO:linopy.io: Writing time: 49.3s
INFO:linopy.solvers:Log file at C:\Users\Mahtab\AppData\Local\Temp\highs.log.
INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 954856 primals, 2119937 duals
Objective: 9.20e+08
Solver model: available
Solver message: optimal

INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-fix-p-lower, Generator-fix-p-upper, Link-fix-p-lower, Link-fix-p-upper, StorageUnit-ext-p_dispatch-lower, StorageUnit-ext-p_dispatch-upper, StorageUnit-ext-p_store-lower, StorageUnit-ext-p_store-upper, StorageUnit-ext-state_of_charge-lower, StorageUnit-ext-state_of_charge-upper, StorageUnit-energy_balance were not assigned to the network.
-24 % Reduction of Renewable Energy installable potential
INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.io:Writing objective.
Writing constraints.: 100%|████████████████████████████████████████████████████████████| 15/15 [00:40<00:00,  2.69s/it]
Writing continuous variables.: 100%|█████████████████████████████████████████████████████| 7/7 [00:07<00:00,  1.08s/it]
INFO:linopy.io: Writing time: 49.45s
INFO:linopy.solvers:Log file at C:\Users\Mahtab\AppData\Local\Temp\highs.log.
INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 954856 primals, 2119937 duals
Objective: 9.20e+08
Solver model: available
Solver message: optimal

INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-fix-p-lower, Generator-fix-p-upper, Link-fix-p-lower, Link-fix-p-upper, StorageUnit-ext-p_dispatch-lower, StorageUnit-ext-p_dispatch-upper, StorageUnit-ext-p_store-lower, StorageUnit-ext-p_store-upper, StorageUnit-ext-state_of_charge-lower, StorageUnit-ext-state_of_charge-upper, StorageUnit-energy_balance were not assigned to the network.
CO2 Shadow price in case of -24 % Reduction of Renewable Energy installable potential:

GlobalConstraint
emission_limit   -7650.545927
Name: mu, dtype: float64

No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.io:Writing objective.
Writing constraints.: 100%|████████████████████████████████████████████████████████████| 15/15 [00:39<00:00,  2.65s/it]
Writing continuous variables.: 100%|█████████████████████████████████████████████████████| 7/7 [00:07<00:00,  1.11s/it]
INFO:linopy.io: Writing time: 49.0s
INFO:linopy.solvers:Log file at C:\Users\Mahtab\AppData\Local\Temp\highs.log.
INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 954856 primals, 2119937 duals
Objective: 2.09e+13
Solver model: available
Solver message: optimal

INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-fix-p-lower, Generator-fix-p-upper, Link-fix-p-lower, Link-fix-p-upper, StorageUnit-ext-p_dispatch-lower, StorageUnit-ext-p_dispatch-upper, StorageUnit-ext-p_store-lower, StorageUnit-ext-p_store-upper, StorageUnit-ext-state_of_charge-lower, StorageUnit-ext-state_of_charge-upper, StorageUnit-energy_balance were not assigned to the network.
1 % Reduction of Renewable Energy installable potential
INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.io:Writing objective.
Writing constraints.: 100%|████████████████████████████████████████████████████████████| 15/15 [00:40<00:00,  2.70s/it]
Writing continuous variables.: 100%|█████████████████████████████████████████████████████| 7/7 [00:07<00:00,  1.14s/it]
INFO:linopy.io: Writing time: 49.95s
INFO:linopy.solvers:Log file at C:\Users\Mahtab\AppData\Local\Temp\highs.log.
INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 954856 primals, 2119937 duals
Objective: 2.09e+13
Solver model: available
Solver message: optimal

INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-fix-p-lower, Generator-fix-p-upper, Link-fix-p-lower, Link-fix-p-upper, StorageUnit-ext-p_dispatch-lower, StorageUnit-ext-p_dispatch-upper, StorageUnit-ext-p_store-lower, StorageUnit-ext-p_store-upper, StorageUnit-ext-state_of_charge-lower, StorageUnit-ext-state_of_charge-upper, StorageUnit-energy_balance were not assigned to the network.
CO2 Shadow price in case of 1 % Reduction of Renewable Energy installable potential:

GlobalConstraint
emission_limit   -4.999777e+06
Name: mu, dtype: float64

No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
In [ ]: